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Abstract 

Accurately  processing  multiple,  time-coincident  signals  presents  a  challenge  to  Elec¬ 
tronic  Warfare  (EW)  receivers,  especially  if  the  signals  are  close  in  frequency  and/or  mis¬ 
matched  in  amplitude.  The  metric  that  quantifies  an  EW  receiver’s  ability  to  measure 
time-coincident  signals  is  the  Instantaneous  Dynamic  Range  (IDR),  defined  for  a  given 
frequency  estimation  accuracy,  a  given  frequency  separation  and  a  given  SNR  as  the  max¬ 
imum  signal  amplitude  ratio  that  can  be  accommodated.  Using  a  two  sinusoid  time-series 
model,  this  thesis  analyzes  IDR  for  ideal  intercept  and  parametric  digital  EW  receivers. 

In  general,  the  number  of  signals  contained  in  the  EW  receiver  measurement  interval 
is  unknown.  Thus,  the  non-parametric  Discrete  Fourier  Transform  (DFT)  is  employed  in 
an  EW  intercept  receiver  with  the  associated  amplitude  dependent  spectral  leakage  which 
limits  IDR.  A  novel  method  to  improve  the  DFT-based  intercept  receiver  IDR  by  com¬ 
pensating  for  the  high  amplitude  signal’s  spectral  leakage  using  computationally  efficient  3 
bin  interpolation  algorithms  is  proposed  and  analyzed.  For  a  desired  frequency  estimation 
accuracy  of  1.5  bins,  the  method  achieves  an  IDR  of  57  dB  with  little  frequency  separation 
dependence  when  the  signals  are  separated  by  more  than  2  bins  with  a  low  amplitude 
signal  SNR  of  10  dB. 

For  situations  where  the  number  of  signals  contained  in  the  measurement  interval  is 
known,  the  IDR  of  an  Rerative  Generalized  Least  Squares  (IGLS)  algorithm-based  para¬ 
metric  receiver  is  analyzed.  A  real  and  complex  signal  IDR  Cramer-Rao  Bound  (IDR-CRB) 
is  derived  for  parametric  receivers  by  extending  results  contained  in  Rife.  For  tight  fre¬ 
quency  estimate  requirements  (these  requirements  depend  on  the  number  of  measurement 
samples),  the  IDR-CRB  yields  achievable  bounds.  For  less  stringent  frequency  estimate 
requirements,  the  IDR-CRB  is  unrealistic  due  to  the  noise  threshold  inherent  to  frequency 
estimation.  Thus,  to  achieve  good  results  when  less  stringent  frequency  estimates  are 
required,  the  author  defines  the  IGLS  algorithm-based  parametric  receiver  IDR  at  the  am¬ 
plitude  ratio  where  the  frequency  estimates  first  achieve  efficiency,  i.e. ,  the  amplitude  ratio 
where  the  overmodelling  condition  first  ceases. 
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ENHANCING  THE  INSTANTANEOUS  DYNAMIC  RANGE  OF  ELECTRONIC 


WARFARE  RECEIVERS  USING  STATISTICAL  SIGNAL  PROCESSING 

I.  Introduction 

The  ability  to  identify  threat  radars  is  of  primary  concern  to  the  warfighter.  When 
flying  missions,  pilots  rely  on  the  Electronic  Warfare  (EW)  system  to  perform  this 
critical  task  by  characterizing  the  threat  radar’s  signal.  If  a  threat  radar  is  nrisidentified  or 
worse,  undetected,  the  consequences  can  be  fatal.  Of  increasing  concern  to  digital  EW  re¬ 
ceiver  designers  is  the  growing  number  of  RF  transmissions;  increasing  the  probability  that 
the  EW  receiver  intercepts  time-coincident  signals.  Because  the  time-coincident  signals  in¬ 
terfere  with  the  receiver’s  ability  to  measure  both  signals  correctly,  special  processing  is 
required  which  is  the  focus  of  this  thesis. 

1.1  EW  Receiver 

EW  receivers  are  unique  from  other  receivers  operating  in  the  RF  region,  e.g.,  radar 
and  communication  receivers.  The  communication  and  radar  receivers  know  the  frequency, 
types  of  modulation,  and  bandwidth  of  the  incoming  signal  [1].  In  an  EW  receiver,  no 
prior  knowledge  of  the  transmission  signal  is  assumed.  In  addition,  in  an  EW  receiver, 
even  the  number  of  intercepted  signals  is  unknown.  Thus,  EW  systems  employ  wideband 
radio  frequency  spectrum  and  signal  analyzers  capable  of  continuous  automatic  real-time 
wideband  search,  detection,  and  analysis  of  signals,  e.g.,  the  Australian  Blue  Owl  System 
[2]- 

An  EW  receiver  measures  certain  signal  parameters  during  each  measurement  interval 
to  enable  the  EW  system  to  identify  the  signals.  Figure  1.1  is  an  illustration  of  the 
parameters,  also  listed  below,  that  are  measured  by  an  EW  receiver. 

•  Frequency 

•  Angle  of  Arrival  (AOA) 
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Figure  1.1:  Parameters  measured  by  EW  receiver  [1], 


•  Pulse  Amplitude  (PA) 

•  Pulse  Width  (PW) 

•  Time  of  Arrival  (TOA) 

•  Polarization  -  EM  polarization  of  pulse  i.e.  vertical,  horizontal,  right  hand  circular, 
left  hand  circular. 

These  parameters  are  used  by  the  EW  system  to  associate  pulses  to  emitters. 

1.1.1  EW  System  Operation.  To  understand  the  purpose  and  requirements  of  an 
EW  receiver,  it  is  also  helpful  to  understand  the  entire  system  operation.  In  this  section, 
a  walk-through  of  the  operations  encountered  by  a  signal  in  the  EW  system  is  discussed. 
The  walk-through  follows  the  block  diagram  of  Figure  1.2  of  a  typical  EW  system  with 
digital  receiver  for  radar  pulse  interception.  The  EW  system  is  analyzed  in  three  parts: 
the  receiver,  the  preprocessor,  and  the  post-processor. 

Following  Fig.  1.2,  the  antenna  intercepts  the  signal  and  propagates  the  signal  to 
the  receiver  RF  amplifier  where  the  information  is  down  converted  to  an  IF  frequency. 
The  down  converted  signals  are  passed  to  the  Analog-to-Digital  Converter  (ADC)  where 
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Simultaneous 


[ 


Digital  Receiver 

Figure  1.2:  Typical  Digital  EW  System  Block  Diagram  [1], 


Table  1.1:  Example  of  a  typical  PDW  format  [1]. 


Parameters 

Range 

No.  of  Bits 

Frequency 

Up  to  32  GHz 

15  (1-MHz  resolution) 

Pulse  Amplitude 

Up  to  128  dB 

7  (1-dB  resolution) 

Pulse  width 

Up  to  204  /as 

12  (0.05  fas  resolution) 

TOA 

Up  to  50  sec 

30  (0.05-/US  resolution) 

AOA 

360  deg 

9  (1-deg  resolution) 

BPSK  signal  flag 

1 

Chirp  signal  flag 

1 

Total  no.  of  bits 

75 

they  are  time  sampled,  i.e.,  a  discrete-time  representation  of  the  signal.  This  discrete-time 
data  is  passed  to  the  Frequency /Spectrum  analysis  module  in  measurement  blocks  where 
spectrum/frequency  estimation  is  performed.  The  frequency/spectrum  analysis  results 
are  passed  to  the  encoder  to  form  the  Pulse  Descriptor  Word  (PDW)  containing  all  of 
the  parameters  for  any  signals  contained  in  the  data.  An  example  of  a  typical  PDW  is 
contained  in  Table  1.1. 

The  preprocessor  processes  the  stream  of  PDW’s  received  from  the  receiver  into 
specific  radar  pulse  trains  through  a  process  called  de-interleaving.  Of  the  five  parameters 
contained  in  a  PDW,  the  three  parameters  used  to  accomplish  de-interleaving  are  the  RF, 
AOA,  and  TOA  difference  between  the  received  pulses  [1].  The  other  two  parameters 
are  unsuitable  for  de-interleaving  because  PA  is  dependent  on  receiving  and  transmitting 
antenna  position  and  PW  is  susceptible  to  multipath  [1] .  Differences  in  AOA  is  by  far  the 
most  stable  parameter  to  use  for  de-interleaving  since  even  aircraft  cannot  quickly  change 
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PW 


Lethal  and  non-lethal 


Frequency 

Figure  1.3:  3-Dimensional  cube  of  parameters  determining  the 

parameter  separating  the  Lethal  threat  from  the  Non-lethal  threat 

[4]- 

their  angular  position  from  pulse  to  pulse  (unfortunately  it  is  also  the  hardest  to  measure 
with  any  accuracy)  [1,3].  Once  the  pulse  trains  are  de-interleaved  into  individual  radar 
pulse  trains,  a  number  of  parameters  can  be  derived  from  the  pulse  train  such  as  antenna 
beamwidth  and  scan  rate  from  successive  amplitude  comparisons,  mode  switching  from 
successive  PW’s,  frequency  of  emitter  pulses  from  multiple  TOA  measurements  referred  to 
as  Pulse  Repetition  Frequency  (PRF),  and  range  from  multiple  AOA  measurements.  This 
information  is  passed  to  the  post-processor  as  an  emitter  report. 

The  post-processor  functions  associate  the  individual  emitter  reports  to  specific  emit¬ 
ters  using  parameters  contained  in  the  emitter  report.  This  parameter  matching  is  analo¬ 
gous  to  determining  where  the  emitter  report  overlaps  with  target  emitters  in  N-dimensions 
as  in  Fig.  1.3  for  a  3-dimensional  cube.  Another  way  to  view  this  process  is  querying  a 
database  with  specific  fields  and  viewing  the  results.  Often,  there  is  overlap  and  multiple 
emitters  match  an  emitter  report.  When  this  occurs  the  emitter  that  is  the  greatest  threat 
is  selected.  When  an  emitter  report  does  not  match  any  entries,  then  an  unknown  emitter 
is  sent.  Once  the  emitter  is  identified,  it  varies  by  platform  on  how  the  information  is 
utilized. 
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1.1.2  EW Receiver  Operation  with  Time-Coincident  Signals.  When  time-coincident 
signals  are  present  in  the  receiver  measurement  interval,  which  is  depicted  in  Fig.  1.4, 
receiver  performance  depends  on  the  spectral/frequency  estimation  function  of  the  re¬ 
ceiver.  If  the  time-coincident  signals  are  characterized  correctly  in  the  spectral  estima¬ 
tion/frequency  estimation  block  of  Figure  1.2,  all  other  systems  will  operate  normally. 
Equi-amplitude  signals  well  separated  in  frequency  do  not  present  much  of  a  problem  for 
the  spectral/frequency  estimation  function.  However,  signals  with  large  amplitude  ratios 
and/or  with  close  frequencies  are  difficult  for  the  receiver  to  measure.  The  receiver  metric 
that  quantifies  the  ability  of  an  EW  receiver  to  measure  time-coincident  signals  is  referred 
to  as  the  Instantaneous  Dynamic  Range  (IDR)  of  an  EW  receiver.  The  standard  IDR 
definition  is 

•  The  standard  IDR  definition  -  The  IDR  is  defined  in  [1]  as  the  power  ratio  of  the 
maximum  and  minimum  simultaneously  received  pulses  that  can  be  properly  encoded 
by  the  receiver  (and  is  similarly  defined  in  [4]). 

Unfortunately,  the  standard  IDR  definition  does  not  reffect  the  dependence  IDR  has  on 
signal  frequency  separation  and  SNR,  which  causes  confusion  when  reporting  results.  Thus, 
the  IDR  definition  employed  in  this  thesis  is 

!  The  thesis  IDR  definition  -  IDR  is  defined  as  the  maximum  signal  amplitude 
ratio  for  a  given  frequency  estimation  accuracy,  a  given  frequency  separation  and 
a  given  Signal-to-Noise  Ratio  (SNR)  [5]. 

1.2  IDR  Analysis  Simplifying  Assumptions 

This  thesis  analyzes  IDR  for  an  ideal  EW  receiver.  The  following  simplifying  as¬ 
sumptions  are  made  for  the  thesis  analysis: 

•  Operation  of  all  devices  prior  to  Spectrum/Frequency  Estimation  block  is  assumed 
nominal  to  include  a  perfect  ADC,  i.e. ,  no  quantization  error. 

•  Signals  are  considered  pure  sinusoidal  tones 

p 

s(n)  =  A*  cos(27r/jn  +  (/>;),  (1.1) 

i—  1 
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Figure  1.4:  EW  Receiver  intercepting  Time- Coincident  Signals. 
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where  s(n)  is  the  received  signal  functional  representation  at  discrete-time  n,  p  is 
the  number  of  tones,  Ai  is  the  ith  signal  (pulse)  amplitude,  ft  is  the  ith  signal  RF 
frequency,  and  (pi  is  the  ith  signal  phase. 

•  Signals  are  assumed  to  fill  the  entire  measurement  period  as  depicted  in  Fig.  1.4. 

•  Thermal  noise  with  the  associated  Additive  White  Gaussian  Noise  (AWGN)  model 
discussed  in  Appendix  A  is  the  only  noise  considered  (no  colored  noise). 

These  same  simplifying  assumptions  are  made  in  [1,5].  Throughout  the  thesis,  signal 
1  is  considered  the  higher  amplitude  signal,  while  signal  2  is  considered  the 
lower  amplitude  signal. 

Under  the  above  simplifying  assumptions,  the  mathematical  measurement  model  is 

x(n)  =  s(n)  +  w(n),  n  =  0, . . , ,N  —  1  (1.2) 

where  the  vector  x=  [x(0),x(l), . . .  ,x(N  —  l)]r  is  the  discrete-time  measurement  (after 
the  ADC  block),  s  is  the  multiple  sinusoidal  tones  defined  in  (1.1),  w  is  AWGN,  and  N  is 
the  number  of  samples.  The  measurement  model  of  (1.2)  is  also  the  standard  model  used 
in  frequency  estimation  in  statistical  signal  processing.  Thus,  there  is  a  vast  amount  of 
literature  discussing  the  analysis  of  (1.2).  Chapter  II  provides  a  brief  review  of  frequency 
estimation  literature  pertinent  to  the  IDR  focus. 

1.3  Problem  Statement 

The  goal  of  this  research  is  to  investigate  the  operation  of  an  ideal  EW  digital  receiver 
when  time-coincident  signals  are  present  to  determine  the  maximum  amplitude  ratio  of  the 
received  signals  at  which  the  receiver  can  still  properly  measure  all  of  the  signals  for  a  given 
frequency  estimation  accuracy,  a  given  frequency  separation  and  a  given  SNR,  referred  to  as 
the  receiver  IDR.  Because,  in  the  ideal  case  considered,  the  measurement  model  is  the  same 
as  in  frequency  estimation;  this  research  applies  statistical  signal  processing  techniques. 
Due  to  the  complexity  of  analyzing  three  or  more  signals,  only  two  time-coincident  signals 
are  considered  in  this  research,  however,  all  techniques  analyzed  (with  the  exception  of  the 
confidence  intervals  established  in  Chapter  V)  can  be  extended  to  three  or  more  signals. 
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EW  receivers  operate  on  real  signals,  thus  real  signals  are  considered  in  the  analysis, 
although  complex  signals  are  sometimes  employed  to  simplify  examples  and  mathematical 
analysis.  The  IDR  definition  employed  throughout  the  thesis  is  defined  in  Section  1.1.2. 

Two  types  of  spectral/frequency  estimation  blocks  are  considered.  The  first  type 
assumes  no  prior  knowledge  of  the  number  of  signals  and  employs  a  non-parametric  sig¬ 
nal  processing  technique  for  spectral/frequency  estimation,  the  Discrete  Fourier  Trans¬ 
form  (DFT).  An  EW  receiver  employing  a  non-parametric  spectral/frequency  estimation 
technique  is  referred  to  as  an  intercept  receiver  in  this  thesis.  The  second  type  of  spec¬ 
tral/frequency  estimation  block  assumes  the  number  of  signals  is  known  and  employs  a 
parametric  based  signal  processing  technique  for  frequency  estimation,  the  novel  Iterative 
Generalized  Least  Squares  (IGLS)  algorithm.  An  EW  receiver  employing  a  parametric 
spectral/frequency  estimation  technique  is  referred  to  as  a  parametric  receiver  in  this  the¬ 
sis. 

As  the  number  of  transmitters  in  the  EW  environment  explodes,  EW  receivers  effec¬ 
tively  processing  time-coincident  signals  is  increasing  in  importance.  Understanding  the 
limitations  and  ways  of  extending  IDR  improves  EW  system  operation,  directly  supporting 
the  war  fighter  in  a  critical  area.  Thus,  this  thesis  has  a  direct  impact  on  USAF  operational 
systems  and  results  can  be  applied  immediately  to  digital  EW  receiver  design  and  software 
updates. 

1.4  Scope 

IDR  is  analyzed  for  the  model  of  (1.2)  using  a  DFT-based  intercept  receiver,  i.e.  a 
non-parametric  frequency  estimator,  and  an  IGLS  algorithm-based  parametric  receiver,  i.e. 
a  parametric  frequency  estimator.  Other  algorithms  are  outside  the  scope  of  this  document. 
The  DFT  is  selected  because  its  universal  applicability  allows  for  use  of  hardware  developed 
for  other  applications  besides  EW.  The  IGLS  algorithm  is  selected  because  it  yields  optimal 
Maximum  Likelihood  frequency  estimates. 
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1.5  Methodology 

The  author  assumes  the  reader  is  familiar  with  Fourier  Transform  (FT)  theory  and 
statistical  signal  processing;  although  Chapter  III  only  requires  FT  theory  and  basic  knowl¬ 
edge  of  statistics.  The  literature  review  in  Chapter  II  is  a  brief  overview  of  applicable 
literature  and  current  IDR  research.  In  lieu  of  reviewing  the  required  mathematical  back¬ 
ground  for  Chapters  III,  IV,  and  V  in  Chapter  II,  each  of  the  Chapters  performs  a  review 
of  the  mathematics  involved  by  (hopefully)  finding  a  common  starting  point  the  reader 
can  follow  in  the  development.  In  Chapter  III,  deterministic  FT  theory  in  the  specific  IDR 
context  is  covered  before  the  theory  is  extended  to  analyze  the  IDR  of  a  novel  multiple 
signal  estimation/detection  technique.  In  Chapter  IV,  the  Cramer-Rao  bound  for  mul¬ 
tiple  equi-amplitude  sinusoidal  signals  originally  derived  by  Rife  in  [6]  is  re-derived,  and 
then  extended  to  IDR  analysis.  In  Chapter  V,  the  novel  frequency  estimation  algorithm, 
Iterative  Generalized  Least  Squares  (IGLS),  originally  developed  by  Dr.  Pachter  and  re¬ 
searched  by  Zahirniak  in  [7]  and  Ingham  in  [8] ,  is  completely  developed  beginning  with  the 
necessary  linear  prediction  background;  concluding  with  an  IGLS  comparison  to  the  IDR 
Cramer-Rao  bound  in  Chapter  IV  and  then  a  modification  to  parametric  IDR  analysis  in 
light  of  the  IDR  Cramer-Rao  bound  comparison  results.  Thus,  original  thesis  results  are 
located  at  the  end  of  Chapters  III,  IV,  and  V  and  are  summarized  in  Chapter  VI.  The 
author  hopes  this  methodology  strikes  a  delicate  balance  between  inundating  the  reader 
with  information  not  pertinent  to  this  research  and  ensuring  the  reader  can  fully  under¬ 
stand  and  interpret  the  research  results.  Finally,  Matlab®  is  employed  where  necessary 
for  simulations  and  analysis. 

1.6  Sponsor 

This  research  is  funded  by  the  Air  Force  Research  Laboratory  Sensors  Directorate 
Radio  Frequency  Analysis  division  Parametrics  branch  (AFRL/SNRP). 

1.7  Summary 

The  research  goal  is  to  investigate  the  performance  of  an  ideal  EW  digital  receiver 
when  time-coincident  signals  are  present  in  the  measurement  using  the  model  of  (1.2) 
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to  determine  the  maximum  amplitude  ratio  of  the  received  signals  at  which  the  receiver 
can  still  properly  measure  all  of  the  signals  for  a  given  frequency  estimation  accuracy,  a 
given  frequency  separation,  and  a  given  SNR,  referred  to  as  the  receiver  IDR.  The  IDR 
is  analyzed  for  a  DFT-based  intercept  receiver  and  an  IGLS  algorithm-based  parametric 
receiver.  This  research  is  accomplished  through  the  use  of  statistical  signal  processing 
techniques  coupled  with  extensive  Matlab®  simulation  and  analysis. 
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II.  Literature  Review 


2. 1  Introduction 

The  literature  pertaining  to  the  thesis  falls  under  three  veins.  Literature  concerning 
IDR,  the  Cramer-Rao  bound,  and  frequency  estimation.  Thus,  all  three  of  these 
topics  are  discussed  below.  This  literature  review  is  intended  to  provide  direction  to 
sources  of  information  to  understand  the  broad  context  of  the  research.  Mathematical 
background  and  literature  pertaining  directly  to  the  methodology  is  covered  in  each  of 
Chapters  III,  IV,  and  V. 

2.2  Instantaneous  Dynamic  Range  Literature  Review 

There  is  little  literature  analyzing  IDR  from  a  statistical  signal  processing  standpoint. 
Most  books  mention  the  issue  and  provide  a  definition,  but  perform  little  analysis  besides 
mentioning  that  IDR  is  frequency  dependent  [1,3].  Part  of  the  problem  is  that  the  standard 
IDR  definition  is  so  broad  that  many  different  interpretations  can  be  inferred.  Some 
engineers  interpret  the  ability  to  measure  the  signal  by  using  a  human  to  interpret  the 
spectrum.  However,  the  EW  system  requires  numerical  frequency  estimates  for  proper 
signal  encoding.  For  this  thesis,  numerical  frequency  estimates  are  required  [9], 
Most  analysis  performed  on  IDR  is  on  a  finished  receiver  with  clarification  on  how  IDR 
is  defined  for  the  tests  seldom,  if  ever,  provided,  i.e.,  whether  the  spectrum  is  interpreted 
or  numerical  estimates  are  obtained,  what  signal  frequency  difference  is  analyzed,  etc.  [9]. 
In  addition  to  clarifying  IDR  for  the  analysis  performed  in  the  thesis,  it  is  hoped  that 
the  tests  performed  in  this  thesis  become  a  standard  for  other  EW  engineers  to  use  when 
reporting  IDR  results. 

2.3  Cramer-Rao  Bound  Literature  Review 

The  standard  optimality  criteria  for  most  estimators  is  Mean  Square  Error  (MSE) 
[10],  defined  as 

MSE  =  T[(6>-0)2],  (2.1) 
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where  9  is  the  parameter  value  to  be  estimated  and  8  is  the  estimated  parameter  value 
(Note  that  RMS  error  is  defined  as  the  square  root  of  mean  square  error).  MSE  is  used 
for  a  variety  of  reasons:  relates  to  a  power  statistic,  penalizes  large  errors  more  than  small 
errors,  and  is  the  error  variance  for  unbiased  estimators  as  shown  below. 

An  illuminating  estimation  result  is  obtained  by  multiplying  out  (2.1)  and  adding 
and  subtracting  £{9}2 


MSE  =  £{82}  -  2 8£{9}  +  92  +  £{9}2  -  £{8}2 
=  var{0}  +  (9  —  £{9})2. 


(2.2) 


Define  bias,  b(9)  as  the  following 


b(9)=£[(8-8)\,  (2.3) 

and  (2.2)  becomes  [10] 

MSE  =  var{<9}  +  b{8)2.  (2.4) 

Thus,  the  MSE  is  composed  of  the  estimate  variance,  var{#},  along  with  the  squared 
estimate  bias. 

Because  of  the  bias  term  of  (2.4),  most  estimators  derived  to  minimize  the  MSE 
directly  are  unrealizable  [10].  However,  estimators  derived  to  minimize  the  estimator 
variance  are  relatively  simple  to  derive,  and  if  the  estimator  can  be  made  unbiased,  the 
estimator  will  minimize  the  MSE  [10].  Limiting  the  analysis  to  unbiased  estimators  also 
allows  comparison  to  the  Cramer- Rao  Bound  (CRB). 

The  CRB  is  a  lower  bound  on  the  error  covariance  matrix  for  any  unbiased  estimator 
of  parameter  6  [11] .  The  CRB  is  standard  fare  in  most  books  on  statistical  signal  processing 
[10-12].  The  CRB  for  multiple  sinusoids  in  white  noise  is  derived  in  [6].  In  [5],  a  Crarner- 
Rao  bound  is  derived  for  complex  signals  IDR,  when  IDR  is  defined  using  the  thesis 
definition,  and  is  referred  to  in  the  thesis  as  the  complex  Instantaneous  Dynamic  Range 
Cramer- Rao  Bound  (IDR-CRB).  The  algorithm  for  the  complex  IDR-CRB  is  simplified  in 
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Chapter  IV  and  extended  to  real  signals.  Note  that  algorithms  that  achieve  the  CRB  are 
referred  to  as  efficient  [11]. 

Unfortunately,  comparison  to  the  IDR-CRB  is  generally  not  valid  for  the  intercept 
receiver  discussed  in  Chapter  III  because  the  DFT  frequency  estimates  are  dominated  by 
frequency  quantization  and  bias  when  multiple  signals  are  present.  Thus,  the  CRB  applies 
mainly  to  the  parametric  receiver  of  Chapter  V,  since  the  estimates  provided  by  the  IGLS 
algorithm  are  unbiased  above  threshold  (the  threshold  effect  is  discussed  in  Chapter  V). 

2-4  Frequency  Estimation  Literature  Review 

Frequency  estimation  is  a  rich  and  varied  subject.  Frequency  estimators  estimate 
the  amplitude,  phase  and  frequency  parameters  of  sinusoidal  signals.  These  parameters 
are  collectively  referred  to  as  the  parameter  vector  9,  where  9  =  [A\  j\  (^1 . . .  Ap  fp  <f>p]T . 
Frequency  estimators  are  divided  into  two  types: 

•  Non-par ametric  Frequency  Estimators:  Frequency  estimators  that  do  not  assume 
any  prior  knowledge  concerning  the  data.  Thus,  non-parametric  techniques  must 
also  determine  the  number  of  signals  present  in  the  measurement  interval 

•  Parametric  Frequency  Estimators:  Methods  that  exploit  the  data  consists  of  the  sum 
of  sinusoids.  Number  of  signals,  referred  to  as  model  order,  is  assumed  known. 

Both  types  of  frequency  estimators  estimate  the  amplitude,  phase  and  frequency  parame¬ 
ters  sinusoidal  signals.  Because  parametric  frequency  estimators  employ  knowledge  of  the 
signal  structure,  the  estimates  provided  are  much  more  accurate  than  the  non-parametric 
frequency  estimators.  However,  if  the  parametric  frequency  estimator  model  order  is  wrong, 
results  are  much  worse  than  non-parametric  frequency  estimators. 

The  DFT  is  the  non-parametric  frequency  estimator  employed  by  the  intercept  re¬ 
ceiver  in  the  thesis  (and  most  digital  EW  receivers).  For  a  single  complex  sinusoid  in  white 
noise,  the  zero-padded  Periodogram  (discussed  in  Appendix  B)  peak  location  is  the  Maxi¬ 
mum  Likelihood  (ML)  estimate  [10,11].  For  one  or  more  real  sinusoids  or  multiple  complex 
sinusoids,  the  location  of  resolvable  Periodogram  peaks  provide  biased  frequency  estimates 
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as  long  as  the  frequencies  are  separated  greater  than  the  Fourier  resolution  (Fourier  reso¬ 
lution  is  defined  in  Chapter  III). 

There  are  many  parametric  frequency  estimation  techniques.  To  limit  the  scope  of  the 
discussion,  only  techniques  that  yield  optimal  ML  frequency  estimates  are  mentioned,  since 
these  are  the  techniques  that  achieve  the  IDR-CRB  above  threshold  (the  interested  reader 
can  refer  to  [13]  for  a  discussion  of  many  sub-optimal  frequency  estimation  techniques). 
The  four  algorithms  found  in  the  literature  that  ML  frequency  estimates  are  listed  below. 

•  The  direct  ML  frequency  estimator  discussed  below. 

•  The  IGLS  algorithm  discussed  in  Chapter  V. 

•  The  IQML  algorithm,  which  is  closely  related  to  IGLS,  discussed  in  [14]. 

•  The  Mean  Likelihood  Frequency  Estimation  algorithm,  discussed  in  [15]. 

2-4.1  Direct  ML  Estimator  for  Complex  Sinusoids  in  AWGN.  The  direct  ML 
frequency  estimate  is  a  complicated  function  of  frequencies  with  many  local  minimum. 
Consider  the  case  of  two  complex  sinusoidal  signals  in  white  noise  (which  is  a  simpler  case 
than  real  signals) 

xc{n)  =  A1ej^hn+^  +  A2ej^hn+*2)  +  w{n),  n  =  0...N-l  (2.5) 

where  the  c  subscript  on  xc  denotes  the  complex  sinusoids.  Form  the  complex  amplitude 
vector  Ac  where  Ac{i)  =  At  exp (jfa).  Now,  let  e*  =  [1,  exp(j2irfi), . . . ,  exp(j2nfi(N  -  1)  )  ]T 
and  the  N  x  2  matrix  E  =  [ei  e2].  The  likelihood  function  of  xc  -  referred  to  as  like¬ 
lihood  because  xc  is  known  and  the  best  estimate  of  the  unknown  parameter  vector 
6  =  [A  1  f\  (f)\ . . .  Ap  fp  <fp\  is  the  6  that  makes  xc  most  likely  to  occur  -  is  [13] 

/e(xc)  =  (7rcr2)_7Vexp  (xc  -  AcE)H(xc  -  ACE)|  .  (2.6) 

If  the  frequencies  are  known,  the  amplitudes  and  phases  ML  estimate  is  a  simple  least 
squares  estimate  [13] 

Ac  =  (E^E)-1E^x.  (2.7) 
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However,  if  the  frequencies  are  also  unknown,  maximizing  the  following  objective  function 
is  required  [13] 

J(f)  =  xHE(EHE)-1EHx.  (2.8) 

J(f)  is  a  highly  non-linear  function  of  the  unknown  frequencies.  This  non-linear  least 
squares  problem  is  computationally  intensive,  and,  if  iterative  techniques  are  employed, 
there  is  no  guarantee  of  achieving  the  global  maximum. 

Figure  2.1  is  a  plot  of  J(f)  with  no  noise  and  signal  parameters:  N=32  (data  record 
length),  [Hi  =  l,_/j  =  0.227,  </>i  =  ^],  [A2  =  l,/2  =  0.207,  ^>2  =  §]■  Note  that  even 
without  the  noise,  the  function  exhibits  many  local  minima  and  maxima  which  complicates 
iterative  maximizing  techniques  -  noise  increases  the  estimation  difficulty  to  the  point 
that  directly  maximizing  J(f)  is  computationally  prohibitive  for  real  time  applications 
[13].  Thus,  most  practical  multiple  signal  frequency  estimation  algorithms  exploit  the 
relationship  between  the  linear  prediction  coefficients  and  frequencies,  which  is  the  basis  of 
the  IGLS  frequency  estimation  algorithm  discussed  in  Chapter  V  and  the  IQML  frequency 
estimation  algorithm  discussed  in  [14]. 

2.5  Conclusion 

This  Literature  review  is  a  brief,  top-level  overview  of  the  issues  involved  in  IDR 
analysis.  IDR,  although  an  important  receiver  concept,  is  normally  only  reported  with 
an  associated  receiver  with  no  discussion  of  how  the  results  are  obtained.  The  CRB  is 
normally  covered  in  a  typical  class  on  statistical  signal  processing,  with  the  CRB  for  mul¬ 
tiple  sinusoids  in  AWGN  derived  by  Rife  in  [6] .  Frequency  estimation,  a  large  and  heavily 
researched  area,  is  discussed  with  emphasis  on  multiple  frequency  estimation  algorithms 
that  achieve  ML  results.  Finally,  in  Chapters  III,  IV,  V,  the  necessary  mathematical 
background,  assuming  the  reader  is  familiar  with  statistical  signal  processing  concepts  and 
Fourier  math,  is  discussed  in  order  to  interpret  the  results  at  the  end  of  the  chapter. 
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(a)  J( f) 


(b)  J(f)  cross-section  at  /i  =  0.227 

Figure  2.1:  Direct  ML  estimator  surface  plot  for  two  complex  sinusoids. 
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III.  IDR  Analysis  of  a  DFT-Based  Intercept  Receiver 
3. 1  Introduction 

In  EW  intercept  receivers,  the  number  of  received  signals  and  the  signals’  carrier  frequen¬ 
cies  are  assumed  unknown.  Thus,  the  intercept  receiver  must  also  detect  the  number  of 
signals  present  in  the  measurement  interval  in  addition  to  measuring  the  signals  frequencies. 
Most  new  digital  EW  receivers  employ  the  DFT  for  multiple  signal  detection/estimation.  If 
the  signals  are  detected,  the  frequencies  are  estimated  from  the  location  of  the  DFT  peaks. 
Thus,  in  this  chapter,  the  DFT-Based  intercept  receiver  IDR  is  analyzed.  In  Section  3.2, 
a  no  noise  deterministic  FT  signal  analysis  is  performed  in  the  context  of  IDR.  Then,  in 
Section  3.3  a  DFT  multiple  frequency  detection/estimation  technique  for  the  DFT-based 
intercept  receiver  is  proposed  and  IDR  evaluated. 

3.2  Deterministic  Analysis:  The  FT  of  the  Sum  of  Two  Sinusoids 

3.2.1  Introduction.  In  deterministic  FT  analysis,  spectral  leakage  limits  IDR. 
Thus,  the  FT  for  the  sum  of  two  sinusoids  is  analyzed  in  this  section  to  fully  understand 
spectral  leakage.  This  section  is  organized  as  follows.  In  Section  3.2.2,  an  infinite  data 
record  is  considered,  which  is  the  easiest  from  an  analysis  standpoint.  Then,  Section  3.2.3 
evaluates  the  impact  of  finite  records  on  the  FT  and  discusses  spectral  leakage.  Data 
record  digitization  is  considered  next  in  Sections  3.2.5  and  3.2.4  including  the  DFT  and 
IDFT;  the  DFT’s  orthogonality  principle  is  also  considered.  While  deterministic  Fourier 
analysis  is  studied  extensively  in  the  literature,  the  subject  matter  covered  in  this  section 
in  the  specific  IDR  context  yields  important  insights  into  the  analysis  of  the  DFT-based 
intercept  receiver’s  IDR.  During  this  sections  analysis,  a  special  type  of  IDR  is  considered, 
called  no  noise  IDR,  which  is  defined  as 

!  no  noise  IDR  -  Maximum  amplitude  ratio  where  the  no  noise  FT  technique 
considered,  i.e.,  FT  or  DFT,  still  exhibits  resolvable  peaks  for  multiple  frequency 
estimates.  This  is  an  analysis  tool  employed  to  aid  the  discussion  on  spectral 
leakage. 
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3.2.2  Infinite  Data  Record.  General  Fourier  theory  is  developed  here  to  motivate 
the  infinite  time  FT  of  two  sinusoids.  If  the  measurement  interval  of  a  signal,  i.e. ,  the  data 
record,  approaches  infinity,  the  signal  can  be  represented  using  the  FT  as  [16] 

/OO 

X(f)e^df  (3.1) 

-oo 

referred  to  as  the  Inverse  Fourier  Transform  (IFT)  of  X(f).  The  FT,  X(f),  represents  the 
frequency  content  of  x(t)  [16] 

/OO 

x{t)e~j2nftdt.  (3.2) 

-OO 

Thus,  x(t)  and  X(f)  are  both  continuous-time  and  continuous-frequency  representations 
of  the  signal  in  the  respective  time  and  frequency  domains. 

The  infinite  FT  of  the  sum  of  two  sinusoids  is  the  simplest  to  analyze  from  an  IDR 
standpoint.  Let  the  infinite  signal,  xj(t),  of  interest  consist  of  the  sum  of  two  carriers, 

X]{t)  =  Aicos(2irfit  +  <fi)  +  A2cos(2tt  f-2t  +  <fo).  (3.3) 


The  FT  of  xj(t)  is 

Aie^1  A^e~^1  Aoe^2  Aoe~^2 

Xi(f)  =  ^^S(f-f1) - l^-S(f  +  f1)  +  -^-S(f-f2)--^-S(f  +  f2).  (3.4) 

Equation  (3.4)  is  the  infinite  time  (ideal)  FT  of  the  sum  of  two  sinusoids  without  noise. 
Note  that  in  the  noiseless,  infinite  time  frequency  representation  of  two  sinusoids,  the 
frequency  resolution  is  theoretically  infinite.  As  far  as  IDR  performance  is  concerned,  the 
frequencies  can  be  detected/estimated  perfectly  regardless  of  frequency  spacing.  Figure  3.1 
illustrates  the  infinite  resolution  for  frequency  detection/estimation  using  the  magnitude 
of  the  FT.  Real  world  finite  signal  records  do  not  yield  infinite  resolution. 

3.2.3  Finite  Data  Record.  For  finite  data  record  analysis,  some  background 
tools  must  be  developed.  To  analyze  the  finite  data  records,  the  rectangular  window  is 
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Figure  3.1:  Infinite  measurement  two  sinusoid  PSD. 


rect  (  —  |  =  < 

T 


introduced 

t\  fl, 

0,  else 

where  r  is  the  measurement  time/period1.  The  FT  of  rect(-)  is 

t\\  cos(7 r/r) 


F  <J  rect  (  -  )  }  = 
r , 


7T/ 


=  rsinc(/r). 


(3.5) 


(3.6) 


Since  the  signals  of  interest  are  always  causal,  the  rectangular  window  is  shifted  in  the 
time  domain  introducing  a  phase  shift  term  in  the  FT  of  the  rectangle  window 


F 


rect  ^  |  =  e  ■77r^rrsinc(/r). 


(3.7) 


The  following  general  convolution  property  of  FT’s  also  aids  the  mathematical  analysis  of 
the  time-limited  signal  [16] 


F{v(t)z(t)}  =  V(f)*Z(f). 


(3.8) 


typical  measurement  time  for  a  digital  EW  receiver  is  100  nanoseconds  [1]. 
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The  above  FT  tools  are  applied  to  the  Fourier  analysis  of  a  signal  consisting  of  the  sum  of 
two  sinusoids. 

The  signal  of  (3.3)  is  multiplied  by  the  time-shifted  rectangular  time  window  to 
model  a  finite  data  record.  The  Fourier  Transform  of  the  finite  data  record  is 

F  (x(t)}  =  F  |x/(f)rect  ^  j>  =  Xj(f )  *  e-J7r^rrsinc(/r),  (3.9) 

using  the  convolution  property  of  (3.8).  Using  the  integral  sifting  property  of  the  5  function 
and  linearity,  the  FT  of  two  finite-time  sinusoids  is 


^(/)=^/(/)«-i7r/Trsinc(/r) 

=  ZLdle_j(^l  +  -2+’r(/+/l)'r)sinc(r(/-|-/1))  — 
+  ld2e-j(^2+7+’r(/+/2)T)sinc(T(/+/2))_ 


Idle— J(— ^l  +  f+7r(/-/l)1')sinc(T(/— /i)) 
ld2e-j(-^2+f+^(/-/2)^)sinc(T(/_/2)) 


(3.10) 


The  sine  function  side-lobe  structure  (called  spectral  leakage  in  academia  or  the  splatter 
effect  by  analog  microwave  receiver  designers)  is  exclusively  a  result  of  the  finite  data 
record  and  limits  the  detection/estimation  of  a  small  amplitude  sinusoid  in  the  presence 
of  a  high  amplitude  sinusoid;  even  in  the  absence  of  measurement  noise.  Side-lobes  can  be 
reduced  using  a  technique  called  windowing.  Windows  are  applied  in  the  same  manner  as 
the  rectangular  function,  and  gradually  reduce  the  transition  of  the  signal  value  to  zero; 
decreasing  the  side-lobe  magnitude  while  increasing  the  main  beam  size  -  Reference  [17]  is 
an  excellent  resource  for  windowing2.  Another  important  limitation  of  the  sine  structure 
besides  sidelobes:  The  signals  cannot  be  detected/resolved  (only  one  peak  will  occur)  if 
the  sinusoids  are  closer  in  frequency  than  ^ ,  referred  to  as  the  Fourier  resolution. 

Figure  3.2  contains  a  plot3  of  the  sum  of  two  sinusoids  FT  with  the  following  pa¬ 
rameters:  [A\  =  1 ,  fi  =  0.32815,  0i  =  |],  [A'2  =  0.01, /2  =  0.1875,02  =  ip],  r  =  32s.  The 
amplitude  ratio  between  the  two  signals  is  40  dB.  In  Fig.  3.2(a),  the  side-lobes  from  the 
higher  magnitude  sinusoid  masks  the  presence  of  the  lower  amplitude  sinusoid.  Thus,  the 

2Note  the  author  has  successfully  applied  windowing  techniques  to  Space-Time  Adaptive  Processing 
STAP)  for  radar  [18]. 

3This  plot  was  generated  using  Matlab’s  DFT  function  and  zero  padding  to  approximate  a  continuous 
FT,  thus  some  minor  aliasing  occurs.  However,  the  general  shape  of  the  continuous  FT  is  the  same. 
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(a)  No  Window 


(b)  Hamming  Window 


Figure  3.2:  Two  Sinusoids  Non-windowed  and  Windowed  FT  Magnitude. 


no  noise  non-windowed  FT  IDR  for  this  particular  frequency  difference  is  less  than  40  dB. 
In  Figure  3.2(b),  the  side-lobes  are  reduced  by  windowing,  and  the  presence  of  the  lower 
amplitude  sinusoid  is  easily  discerned.  The  no  noise  Hamming  windowed  FT  IDR  at  this 
particular  frequency  difference  is  greater  than  40  dB. 

Figure  3.3  contains  a  plot  of  the  sum  of  two  sinusoids  FT  with  the  following  param¬ 
eters:  [Hi  =  1,/i  =  0.32815,0!  =  f],  [A'2  =  1  ,/2  =  0.1875,02  =  ^] ,  r  =  32s.  These 
frequencies  are  within  the  Fourier  resolution,  thus  the  two  frequencies  cannot  be  re¬ 
solved  using  conventional  FT  methods  (windowing  exacerbates  the  situation  because  of 
the  widening  of  the  main  beam).  The  no  noise  FT  IDR  for  these  two  signals  is  zero;  they 
are  not  resolved. 


3.2.4  Sampling.  Since  the  EW  receiver  samples  the  signal,  the  sampling  effects 
on  the  above  developments  are  now  analyzed;  however,  quantization  effects  are  ignored. 
Assuming  the  signal  is  sampled  uniformly,  the  following  sampling  function  is  introduced 
to  model  the  sampling  [19] 

OO 

g(t)  =  Y  S(t-nT),  (3.11) 

n=— 00 
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Figure  3.3:  Two  Sinusoids  within  Fourier  Resolution  FT  Magnitude. 

where  T  is  the  sampling  period.  Equation  (3.11)  is  periodic  and  can  be  represented  in 
terms  of  a  Fourier  Series  [19], 


g(t)=  £  Cn 


Jn2nf0t 


(3.12) 


where  fa  is  the  sampling  frequency,  f0  =  b,  and  Cn  is  defined  as 


Cn  =  f  S(t)e-^otdt 


(3.13) 


Recall  that  x(t)  is  the  signal  of  interest  and  let  xs(t)  represent  the  sampled  version  of  the 
signal, 


xs(t)  =  g(t)x(t) 


E 


(3-14) 
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Taking  the  FT  of  (3.14)  and  reversing  the  order  of  integration  and  summation  yields  the 
following  expression 

/OO  00  I 

^2  -ejn2nfotx{t)e-j2nftdt 

n=— 00  T 
1  _°°_  /*oo 

=  -  ^2  /  x(t)e-j2*{f~nfo)tdt  (3.15) 

1  n=- 00  ^-°° 

1  00 

=  T  E  XU~nfo). 

72—  —  OO 

Thus,  as  (3.15)  illustrates,  the  sampled  signal’s  FT  is  a  periodic  version  of  the  continuous 
signal  FT  about  the  sampling  frequency.  The  major  developments  for  the  continuous  two 
sinusoid  signal  will  hold  for  the  sampled  sinusoid  signal  as  long  as  the  sampling  rate  is  at 
least  twice  the  highest  frequency  in  the  signal,  to  prevent  aliasing.  However,  since  finite 
signals  have  infinite  frequency  content,  the  infinite  side-lobes  of  the  sine  function  will  cause 
some  distortion  effects  from  aliasing  in  (3.10). 

Figure  3.4  contains  a  plot  of  a  sampled  signal  FT  and  a  continuous  FT  for  a  single 
sinusoid  with  the  following  simulation  parameters:  [Hi  =  1,  f\  =  0.32815,  (j)\  =  |],  N  =  32. 
The  effects  of  aliasing  are  noticeable  on  side-lobes  near  /  =  0.5  and  /  =  —0.5.  The  aliasing 
causes  the  side-lobes  of  the  sampled  signal  FT  to  be  higher  than  the  continuous  FT,  while 
effects  near  the  main  lobe  are  negligible. 

3.2.5  Discrete  Fourier  Transform  of  Two  Sinusoids.  Since  the  intercept  receiver 
is  a  digital  receiver,  the  frequency  domain  is  also  digital.  The  FT  of  a  finite,  uniformly 
sampled  signal  is 

N-l 

X{f)  =  J2  x{nT)e~j2nfnT .  (3.16) 

71=0 

The  signal  is  uniformly  sampled  in  the  frequency  domain  by  letting  /  =  k  = 

0, . ,N  —  1  [19].  Equation  (3.16)  then  becomes 

N-l 

X(k )  =  ^2  x(nT)e~j27TT? ,  (3.17) 

77=0 
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Figure  3.4:  Single  sinusoid  continuous  and  sampled  FT  magnitude. 

where  X(k)  is  the  DFT  of  x(nT )4.  X(k )  is  referred  to  as  the  Discrete  Fourier  Spectrum 
(DFS).  Note  that  if  the  number  of  frequency  samples  is  a  power  of  two  (the  number  of 
frequency  samples  can  be  more  than  the  number  of  time  samples  if  the  signal  is  zero 
padded),  the  number  of  arithmetic  operations  required  to  perform  the  DFT  is  reduced 
dramatically  using  the  Fast  Fourier  Transform  (FFT)  algorithm.  The  inverse  of  the  DFT, 
the  Inverse  Discrete  Fourier  Transform  (IDFT),  is 

N- 1 

F-1  {X(k)}  =  x(n)  =  -  Y,  X(k)e^.  (3.18) 

k= 0 

Through  use  of  the  IDFT  and  DFT,  a  discrete-time  signal  can  be  digitally  represented  in  the 
time  or  frequency  domain.  The  resolution  in  frequency  when  the  DFT  is  used  is  delimited 
by  the  number  of  frequency  sample  points,  i.e.  frequency  quantization  [15].  Increasing  the 
frequency  sampling  by  zero  padding  the  time  sequence  increases  the  frequency  resolution, 
but  at  a  cost  of  increased  computational  complexity.  When  the  sinusoids  are  orthogonal, 

4Further  developments  will  assume  that  the  sampling  rate,  T,  is  equal  to  one  unless  stated  otherwise 
since  all  results  can  be  scaled  for  the  appropriate  sampling  rate  (This  is  a  standard  practice  in  most  of  the 
literature  on  Signal  Processing). 
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(a)  No  Zero  Padding  (b)  Zero  Padding  to  64  Bins 

Figure  3.5:  Two  Sinusoids  Zero  Padded  and  Not  Zero  Padded  DFT  Magnitude. 

i.e.  their  frequencies  are  exactly  sample  points  of  the  DFT,  the  side-lobes  do  not  mask  the 
signals,  and  the  no  noise  DFT  IDR  is  infinite. 

Figure  3.5  contains  a  plot  of  the  sum  of  two  sinusoids  zero  padded  and  non-zero 
padded  DFT  with  the  following  parameters:  [A\  =  l,/i  =  0.32815,  0i  =  |],  [A2  =  1,  f-2  = 
0.1875,  fa  =  xl’  N  =  32.  The  outline  of  the  continuous  FT  (generated  using  a  4096- 
point  zero-padded  DFT)  is  also  plotted  to  illustrate  where  the  samples  are  located.  Fig. 
3.5(a)  shows  the  frequency  sample  points  are  referred  to  as  bins,  numbered  0  to  31  starting 
with  the  positive  frequencies.  A  common  way  to  refer  to  frequency  in  the  DFT  is  the  bin 
position.  Thus,  f\  =  0.32815  is  bin  position  10.5,  f\  =  10.5/M.  The  plot  in  Fig.  3.5(b) 
illustrates  the  point  that  zero  padding  increases  the  number  of  frequency  samples. 

Figure  3.6  contains  of  plot  of  the  sum  of  two  sinusoids  for  orthogonal  and  non- 
orthogonal  sinusoids  over  the  observation  interval  for  the  simulation  parameters:  [A  \  = 
200,  f\  =  Bin  10,  </>i  =  |],  [A2  =  2, /2  =  Bin  6,02  =  xl>  ^  =  32.  For  the  orthogonal 
sinusoids  plotted  in  Fig.  3.6(a),  the  no  noise  DFT  IDR  is  infinite  in  because  the  frequency 
samples  lie  along  the  zeroes  of  the  other  signals  sine  functions.  For  Fig.  3.6(b),  the 
frequency  of  f\  is  changed  to  Bin  10.5.  The  spectral  leakage  of  signal  1  is  maximum  in 
this  case,  and  signal  2  cannot  be  resolved.  Thus  in  this  situation,  the  no  noise  DFT  IDR 
is  less  than  40  dB. 
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(a)  Orthogonal  (b)  Not  Orthogonal 


Figure  3.6:  Orthogonal  and  Not  Orthogonal  Sinusoids  DFT  magnitude. 

3.2.6  Conclusion.  The  above  discussion  shows  that  Fourier  techniques  are  quite 
robust  for  the  detection  and  estimation  sinusoids  at  least,  so  far,  in  the  absence  of  noise. 
The  two  areas  where  the  no  noise  FT  needs  compensation  to  increase  no  noise  IDR  are 

1.  Sinusoids  within  the  Fourier  Resolution  Limit.  Other  frequency  estimation 
and  detection  methods  must  be  employed  when  within  the  Fourier  resolution,  re¬ 
gardless  of  amplitude. 

2.  Detection  interference  from  the  side  lobes  of  a  stronger  signal.  Windowing 
is  one  popular  method  to  increase  no  noise  IDR. 

The  next  section  analyzes  a  DFT-based  EW  intercept  intercept  receiver  for  its  IDR  in  the 
presence  of  sidelobes  and  noise. 

3.3  DFT-Based  Intercept  Receiver  Frequency  Detection/ Estimation  in  Noise 
Now,  return  to  the  measurement  model  of  (1.2) 

x(t)  =  Ai  cos(27t/i  +  <fi)  +  A2  cos(27t fin  +  <f/)  +  w(n).  (3.19) 

In  the  following  analysis,  signal  1  is  considered  the  higher  amplitude  signal,  while  signal 
2  is  considered  the  lower  amplitude  signal.  Determining  the  number  of  sinusoids  present 
for  signals  buried  in  white  noise  is  a  difficult  problem,  especially  for  high  amplitude  ratio 
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signals.  The  most  accurate  methods  employ  statistical  tests  that  attempt  to  find  the 
hypothesis  with  the  least  error  energy,  such  as  information  theoretic  criteria  approaches 
[13].  Such  methods  are  preferred  because  they  avoid  the  use  of  thresholds  and  have  higher 
resolution  than  the  Fourier  resolution.  However,  the  methods  are  computationally  intensive 
and  not  practical  for  Electronic  Warfare  applications. 

The  most  robust  method  for  an  EW  intercept  receiver  is  to  use  either  the  DFS  or 
Periodogram  -  described  in  Appendix  B.  Assuming  the  noise  variance  is  known,  a  noise 
threshold  can  be  set  for  detection  by  determining  an  acceptable  probability  of  false  alarm 
Pfa,  i.e.  probability  that  the  noise  will  exceed  the  threshold  at  any  frequency  bin  for  a 
single  measurement  interval,  and  setting  the  threshold  accordingly5.  If  any  bin  exceeds  the 
threshold,  a  detection  is  declared,  and  the  maximum  bin  of  the  measurement  is  a  coarse 
(quantized)  frequency  estimate  (if  noise  exceeds  the  threshold,  it  is  a  false  alarm). 

Figure  3.7  contains  a  demonstration  of  the  above  concepts.  The  noise  alone  exceeds 
the  bin  threshold  for  Fig.  3.7(a),  thus  a  false  alarm  would  be  declared  for  this  measurement 
interval.  If  multiple  bins  exceed  the  threshold  for  a  given  measurement,  only  one  false  alarm 
is  declared.  In  Figure  3.7(b),  a  single  sinusoid  has  exceeded  the  threshold  and  a  detection 
is  declared  with  simulation  parameters6:  [A\  =  l,/i  =  10.5  bins,^>i  =  -y],  SNR  =  20 
(dB),  N  =  32.  A  coarse  frequency  estimate  for  the  sinusoid  is  the  bin  location  of  the 
maximum  peak,  as  shown  in  the  figure.  Note  that  for  the  assumed  single  sinusoid  case, 
the  side-lobes  exceeding  the  threshold  are  ignored. 

If  more  than  one  signal  could  be  in  the  measurement  record,  as  is  the  case  in  EW, 
the  detection  problem  becomes  more  complex  because  of  the  spectral  leakage.  Fig.  3.8 
contains  a  plot  of  the  DFT  with  simulation  parameters:  [A±  =  1,  fi  =  10.5  bins,^>i  =  |], 
[Ar2  =  0.2, /2  =  14.2  bins,(/>2  =  ^],  SNR,  =  20  (dB)  referenced  to  signal  1,  N  =  32. 
Note  that  the  side-lobe  bin  magnitude  pointed  to  in  the  figure  is  above  the  signal  2  bin 
magnitude.  EW  receivers  require  that  the  presence  of  the  weaker  signal  must  be  recognized 
by  automation;  the  spectral  leakage  exceeding  the  threshold  near  the  main  lobe  must  be 

sTypical  Pfa  for  an  EW  receiver  is  10-11. 

6 Since  real  signals  are  considered,  from  here  on  out  in  this  development  only  bins  1  through  4f  (sample 
number  is  assumed  even)  will  be  considered;  the  negative  frequency  images  will  be  ignored  unless  stated 
otherwise. 
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(a)  Noise  (b)  Noise  +  1  Sinusoid 


Figure  3.7:  Threshold  Demonstration  with  Noise  and  Noise+1  Sinusoid. 

disregarded,  and  the  presence  of  the  weaker  signal  must  be  recognized.  This  is  a  difficult 
pattern  recognition  problem  to  implement  in  a  machine  [1].  Recall  that  the  frequency 
estimate  must  be  numerical.  There  is  no  human  evaluating  a  DFS  display. 

A  method  employed  to  avoid  the  pattern  recognition  problem  by  the  AFRL  Labo¬ 
ratory  exploits  the  prior  knowledge  of  the  spectral  leakage  shape  to  subtract  the  spectral 
leakage  [9].  Thus,  a  method,  referred  to  here  as  the  Spectral  Leakage  Reduction  (SLR) 
method,  with  the  same  general  concept  is  proposed  and  analyzed  here.  The  method  entails 
the  following  steps 

1.  Find  the  peak  of  the  DFS  and  ensure  it  is  above  threshold. 

2.  Estimate  the  high  amplitude  signal’s  frequency,  amplitude,  and  phase. 

3.  Subtract  the  spectral  leakage  from  the  DFS  (this  subtraction  would  be  accomplished 
via  a  lookup  table  with  the  actual  system;  for  simulations,  the  DFS  is  calculated). 

4.  Find  the  peak  of  the  subtracted  DFS  and  check  if  it  is  above  threshold. 

5.  Declare  a  detection  if  the  peak  is  above  threshold. 

If  the  sinusoidal  model  is  correct,  the  above  method  hinges  on  the  accuracy  of  the  frequency, 
phase,  and  amplitude  estimates  for  accurate  spectral  leakage  estimation.  This  method 
is  similar  to  a  method  proposed  in  [20],  although  [20]  is  using  the  method  to  reduce 
bias  for  extremely  accurate  interpolation  multiple  frequency  estimation  (the  approximate 
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Figure  3.8:  Two  Sinusoid  +  Noise  DFT  Magnitude. 

locations  of  the  peaks  are  assumed  known).  Whereas,  the  SLR  method  is  reducing  the 
spectral  leakage  in  order  to  detect  low  amplitude  sinusoids  in  the  presence  of  high  amplitude 
sinusoids. 

Figure  3.9  contains  a  plot  of  the  signal  DFT  of  Figure  3.8  with  the  spectral  leakage 
of  signal  1  compensated  for  by  subtracting  the  DFT  magnitude  of  signal  1  from  the  DFT 
magnitude  of  signal  1  and  2  +  noise 


XsLR(k)  =  I  S2(k)  +  S^k)  +  W{k)  -  Si  (A;)  I 
=  \S2(k)  +  W(k)  +  eSLR(k)\ 


(3.20) 


where  XsLR^k)  is  the  SLR  method  DFS  (with  the  spectral  leakage  from  signal  one  com¬ 
pensated  for),  S\(k)  is  the  estimated  DFS  of  signal  1,  S\(k)  is  the  estimated  DFS  of  signal 
1,  W(k)  is  the  DFS  of  the  noise,  S2(k)  is  the  estimated  DFS  of  signal  2  and  e.SLR  is  the 
error  between  the  estimated  signal  1  DFS  and  the  actual  DFS  of  signal  1.  In  Fig.  3.9(a), 
the  first  signal  frequency  is  known  exactly  and,  as  expected,  the  compensation  removes 
the  side-lobes  to  below  threshold,  allowing  the  second  signal  to  be  detected  via  a  threshold 
test,  i.e.  esLR  is  very  low.  In  Fig.  3.9(b),  the  frequency  used  to  estimate  the  side-lobes  is 
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(a)  Compensated  exact  frequency.  (b)  Compensated  0.25  bins  difference  frequency 

Figure  3.9:  Sine  Lobe  Compensation  Demonstration. 


at  10.75  bins  and,  as  a  result  of  the  poor  frequency  estimate,  the  side-lobes  are  not  ade¬ 
quately  compensated  for,  i.e.  a  larger  esLR ■  Thus,  good  frequency  estimates  are  required 
for  the  SLR  method. 

3.3.1  Frequency  Interpolation  Methods.  DFS  peak  interpolation  algorithms  are 
developed  in  this  section  in  order  to  estimate  the  frequency  for  the  SLR  method  accurately. 
The  algorithms  analyzed  are  computationally  efficient  and  only  rely  on  the  three  DFS 
points  closest  to  the  actual  sinusoids  frequency  (note  that  these  three  points  contain  >  85% 
of  the  signal  energy  [20]),  thus  they  are  suitable  for  EW  applications.  All  algorithms  are 
developed  for  the  one  complex  sinusoidal  signal  condition,  or  cisoid,  because  there  is  no 
spectral  leakage  bias  to  contend  with  (for  one  real  signal,  the  negative  frequency  image  has 
spectral  leakage  in  the  positive  frequency  image) .  Simulations  are  performed  to  gauge  the 
effect  of  the  presence  of  one  and  two  real  sinusoids  on  algorithm  performance. 

3.3. 1.1  Modulus  Peak  Position  Interpolation.  The  first  DFS  interpolation 
algorithm,  called  here  the  Modulus  Peak  Position  (MPP)  interpolator,  is  discussed  in  [1,20]. 
The  MPP  estimates  the  points  by  using  the  largest  DFS  absolute  value  peak  and  the  largest 
absolute  value  of  the  peak’s  two  neighbors  (the  one  also  lying  on  the  main  lobe).  It  then 
interpolates  the  value  in  between  them.  The  mathematics  of  this  interpolation  follows. 
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First,  find  the  maximum  DFS  sample  point  magnitude,  max(|X(A;)|)  equals  |X(p)|. 
Next,  find  the  neighboring  peak  with  the  highest  absolute  value  amplitude  and  set  the 
variable  a  according  to  the  position,  that  is 


a  =  < 


1,  if  \x(p  +  1)  |  >  \x(p  —1)|; 
—  1,  otherwise. 


(3.21) 


Now,  assume  the  DFS  value  is  X(^)  where  A  is  the  distance  from  the  true  signal  frequency 
in  fractions  of  a  bin  jj{p  +  A)  =  /,  where  A  <  |0.5|.  Thus,  the  value  is 


.  A\  .  f  A\  sin(7rA) 

x  (  T7  )  =  Npsmc  (  —  )  =  P - 


N 


N 


(3.22) 


where  P  is  the  continuous  signal  complex  amplitude.  The  highest  magnitude  neighbor  DFS 
value  is 


(3.23) 


N 


Now,  A  is  interpolated  as  [1,20] 


A  =  a 


|  X(m  +  a)  | 


|A(m)|  +  | X(m  +  a) | ' 


(3.24) 


The  corresponding  frequency  estimate  is 


/  =  (jp  +  A)  —  . 


(3.25) 


where  the  sampling  frequency  is  assumed  one  for  the  above  development. 

Figure  3.10  is  an  illustration  of  the  above  math.  For  larger  A,  as  in  Fig.  3.10(a), 
the  amplitude  difference  is  well  pronounced  between  | X(p  —  1)|  and  \X(p+  1)|  (For  this 
example  a  =  —1).  For  small  A,  as  in  Fig.  3.10(b),  the  amplitude  difference  is  negligible 
and  this  will  cause  errors  when  noise  is  present. 

Small  systematic  errors  are  present  in  all  of  the  interpolation  algorithms  (including 
the  above  method)  which  cause  a  slight  bias  in  the  estimate  [20].  When  real  signals  are 
considered,  the  negative  frequency  image  side-lobes  introduce  larger  bias  to  the  interpolator 
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(a)  Large  A.  (b)  Small  A. 


Figure  3.10:  MPP  concept  illustration. 

frequency  estimates.  Figure  3.11  contains  a  MPP  no  noise  frequency  estimate  deterministic 
bias  analysis  for  a  single  real  sinusoid  and  complex  cisoid  with  the  following  parameters: 
[A i  =  1,  (f>i  =  0].  In  the  figure,  the  deterministic  bias  is  calculated  between  the  interpolated 
frequency,  / ^ ,  and  the  actual  frequency,  /a,  for  the  bin  positions  indicated 

bias  =  |  /a  —  /a  I  (3-26) 

at  0.01  intervals  of  A  and  plotted  as  — 10  log10(bias)  (which  means  good  performance  is 
plotted  above  bad  performance) .  The  negative  frequency  side- lobe  aliasing  into  the  positive 
frequencies  causes  an  increase  in  bias  of  approx  30  dB  at  some  points  over  the  complex 
cisoid  case.  This  real  signal  bias  value  varies  according  to  the  bin  position.  This  real  signal 
bias  would  be  reduced  by  using  a  window  based  interpolator  because  of  the  associated 
spectral  leakage  reduction,  but  windowing  adversely  effects  probability  of  detection  and 
frequency  resolution  of  close  frequency  signals. 

Figure  3.12  contains  a  pseudocolor  plot  of  a  Monte  Carlo  (MC)  MSE  analysis  for 
estimating  a  real  sinusoid’s  frequency  using  the  MPP  interpolator  with  one  real  sinusoid 
and  noise,  where  MSE  is  calculated  as 

1  M 

MSEvar  =  -YJU^-fl?-  (3-27) 

i= 1 
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Figure  3.11:  Deterministic  bias  of  MPP  interpolator  for  real  and  complex  signals. 

where  M  is  the  number  of  Monte  Carlo  trials  and  is  the  ith  frequency  estimate  of  the 
MC  trial.  For  the  intercept  receiver  simulations,  the  phase  is  randomly  distributed  between 
0  and  2ir.  The  plot  is  generated  by  performing  1000  MC  trials  at  each  of  33  evenly  spaced 
A  from  bin  32  and  repeating  for  each  SNR  step  with  simulation  parameters:  [A\  =  1,  <f>\  = 
C/[0,27t]],  N=256.  The  color  bar  indicates  the  — 101og10(MS,£'),  thus  the  higher  the  value, 
the  better  the  estimate.  From  the  plots,  the  MPP  interpolator  performance  significantly 
degrades  for  low  A  values.  This  poor  performance  is  due  to  the  a  sign  being  switched  at 
low  A,  which  introduces  a  larger  error  and  bias.  Also,  the  bias  from  the  real  signal  spectral 
leakage  dominates  the  MSE  above  30  dB  SNR,  which  is  why  the  MSE  values  plateau  at 
this  point. 

3.3. 1.2  Phase  Based  Interpolator.  Because  of  inaccuracy  in  determining  a 
for  low  A,  the  following  phase  based  interpolation  algorithms,  called  the  Phase  Based  In¬ 
terpolator  (PBI),  is  introduced  to  overcome  this  deficiency  [20].  First,  the  two  neighboring 
DFS  bins  phase  X(a)  are  referenced  to  the  largest  bin’s  phase  and  then  the  real  part  of 
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Figure  3.12:  MPP  Interpolator  MSE  versus  A  and  SNR.  Color  bar  indicates  value  of 

—  101og10 (MSE).  Thus  the  higher  the  value  is,  the  better  the  performance. 


this  operation  is  taken  [20] 


V(fi)  =  Re{X(p  +  n)X*(p)}  (3.28) 

where  V(fi)  is  called  the  Phase  Indexed  Variable  (PIV).  The  following  test  is  used  to 
calculate  A  [20] 

iff  V(-l)  -  V(l)  >  0,  A  =  A+,  elseA  =  A_.  (3.29) 

where  A+  and  A_  are  defined  as  [20] 

A  ~r(i)  A  v(-i )  ,,,a) 

+  V(0)  —  V(l)  “  F(0)  —  F(— 1)  1  '  ’ 

As  A  approaches  zero,  V(l)  ~  V(— 1)  and  the  PBI  estimate  exhibits  much  less  error  than 
the  MPP. 

Figure  3.13  contains  a  PBI  no  noise  frequency  estimate  bias  analysis  for  a  single  real 
sinusoid  and  complex  cisoid  with  the  same  setup  as  Fig.  3.11.  The  negative  frequency 
sidelobe  aliasing  into  the  positive  frequencies  causes  an  increase  in  bias  of  approx  20  dB  at 
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Figure  3.13:  Deterministic  bias  of  PBI  interpolator  for  real  and  complex  signals. 

some  points  over  the  complex  cisoid  case.  This  real  signal  bias  value  varies  according  to  the 
bin  position,  which  is  no  different  than  the  MPP.  The  PBI  complex  cisoid  bias  is  slightly 
higher  than  the  MPP  bias.  However,  the  PBI  does  not  experience  the  jagged  increase  in 
bias  at  low  A. 

Figure  3.14  contains  a  pseudocolor  plot  of  a  MC  MSE  analysis  for  estimating  a  real 
sinusoid’s  frequency  using  the  PBI  interpolator  with  one  real  sinusoid  and  noise  with  the 
same  setup  as  Fig.  3.12.  From  comparing  Fig.  3.14  to  Fig.  3.12,  the  PBI  performance 
significantly  improves  for  low  A  values  over  the  MPP,  which  is  expected  because  of  the 
phase  referencing  to  decrease  the  increased  MPP  bias  associated  with  picking  the  wrong 
a  value.  Also,  the  bias  from  the  real  signal  spectral  leakage  dominates  the  MSE  above  30 
dB  SNR,  which  is  why  the  MSE  values  plateau  at  this  point  (except  for  the  spike  locations 
in  3.13  for  bin  32). 

3.3. 1.3  Gamma  Phase  Based  Interpolator.  The  PBI  can  also  be  improved 
upon  by  using  the  following  algorithm,  the  Gamma  Phase  Based  Interpolator  (GPBI). 
When  A  is  small,  V(l)  and  V(-l)  provide  independent  estimates  of  A.  Thus,  some  esti¬ 
mation  gain  over  PBI  is  gained  by  averaging  the  two  [20].  In  [20],  the  following  average  is 
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Figure  3.14:  PBI  Interpolator  MSE  versus  A  and  SNR.  Color  bar  indicates  value  of 

— 10  logw(MSE).  Thus  the  higher  the  value  is,  the  better  the  performance. 

proposed 

n-i)  -  ni) 

'  2V  (0)  +  V(— 1)  +  V(l) 

Using  7,  Reference  [20]  gives  the  A  estimate  as 

A_  \/l  +  872  —  1  . 

47 

Figure  3.15  contains  a  GPBI  no  noise  frequency  estimate  bias  analysis  for  a  single 
real  sinusoid  and  complex  cisoid  with  the  same  setup  as  Fig.  3.11.  The  negative  frequency 
side- lobe  aliasing  into  the  positive  frequencies  bias  is  reduced  from  the  PBI  and  MPP  by 
using  the  gamma  interpolator.  The  GPBI  complex  cisoid  bias  is  also  slightly  higher  than 
the  MPP  cisoid  bias.  Again,  however,  the  GPBI  does  not  experience  the  jagged  increase 
in  bias  at  low  A. 

Figure  3.16  contains  a  pseudocolor  plot  of  a  MC  MSE  analysis  for  estimating  a  real 
sinusoid’s  frequency  using  the  PBI  interpolator  with  one  real  sinusoid  and  noise  with  the 
same  setup  as  Fig.  3.12.  From  the  plots,  the  GPBI  shows  some  improvement  over  the 
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Figure  3.15:  Deterministic  bias  of  GPBI  interpolator  for  real  and  complex  signals. 

PBI  and  marked  improvement  over  the  MPP.  Also,  the  bias  from  the  real  signal  spectral 
leakage  dominates  the  MSE  above  30  dB  SNR,  which  is  why  the  MSE  values  plateau  at 
this  point. 


3.3. 1-4  Two  Sinusoid  Interpolator  Performance.  With  two  sinusoids,  the 
interpolator  estimate  bias  increases.  Figure  3.17  is  a  two  sinusoid  bias  analysis  of  the 
interpolator  algorithms  for  the  following  parameters:  [A\  =  1,  <p\  =  0],  [A2  =  0.8,  ^>2  =  0]. 
Three  bin  positions  are  considered  for  the  first  sinusoid:  32,  64,  and  100.  The  second 
sinusoid  is  at  bin  35.5  for  maximum  leakage.  Sinusoids  within  2  bins  will  be  considered  high 
resolution  for  the  DFT-based  intercept  receiver  and  are  not  considered  in  this  Chapter  [9] . 
For  reference,  the  bias  of  using  |A(p)|  as  an  estimate  is  plotted.  For  the  MPP  estimate 
in  Fig.  3.17(a),  the  low  A  performance  decreases  considerably;  below  using  the  FFT 
bin  |X(p)|  as  the  estimate.  The  PBI  and  GPBI  for  Figures  3.17(b)  and  3.17(c)  both 
outperform  using  the  FFT  bin  as  the  estimate,  with  the  GPBI  performing  the  best  of  the 
three  interpolation  algorithms.  In  terms  of  IDR,  the  2  sinusoid  interpolator  bias  decreases 
as  the  IDR  increases,  which  is  good  in  terms  of  estimating  the  spectral  leakage.  In  other 
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Figure  3.16:  GPBI  Interpolator  MSE  versus  A  and  SNR.  Color  bar  indicates  value  of 

—  10  log i0(MSE).  Thus  the  higher  the  value  is,  the  better  the  performance. 


words,  the  estimate  performance  improves  as  the  IDR  increases,  which  in  turn  means  the 
spectral  leakage  estimate  is  improving,  which  is  a  desirable  characteristic. 

Figure  3.18  contains  a  pseudocolor  plot  of  a  two  sinusoid  MSE  analysis  of  the  in¬ 
terpolator  algorithms  with  the  same  signal  parameters  as  Fig.  3.17  except  the  phase  is 
uniformly  distributed  for  both  signals  between  0  and  2n.  The  colorbar  indicates  the  value 
of  —  101og10(MS'£')  between  /i  and  f\  for  1000  MC  runs  at  the  specified  SNR  and  A  value. 
The  first  sinusoid  is  centered  at  bin  32  with  34  A  evenly  distributed  around  the  bin.  The 
second  sinusoid  is  located  at  bin  35.5  for  maximum  leakage  with  the  amplitude  ratio  of 
the  signals  maintained  the  same  as  Fig.  3.17.  The  close  proximity  of  frequency  and  ampli¬ 
tude  of  the  first  signal  to  the  second  signal  makes  this  a  demanding  test.  As  expected,  the 
MPP  algorithm  exhibits  the  worst  MSE  performance  of  the  three  interpolation  algorithms, 
especially  for  small  A.  The  GPBI  algorithm  exhibits  a  slightly  better  performance  than 
the  PBI.  Above  20  dB  SNR,  the  MSE  of  all  of  the  algorithms  is  dominated  by  the  bias. 
Once  again  though,  the  estimate  performance  improves  as  the  IDR  increases,  which  in  turn 
means  the  spectral  leakage  estimate  is  improving,  which  is  a  desirable  characteristic. 
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Figure  3.18:  Interpolator  2  Sin  interpolator  estimate  f\  MSE  versus  A  and  SNR.  Col- 
orbar  indicates  —  lOlogio(MSE)  value. 
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3-4  SLR  Method  Analysis 


As  discussed  in  the  introduction  to  this  section,  the  following  algorithm,  called  the 
SLR  Method,  is  proposed  to  increase  weak  signal  detection/estimation  performance  using 
the  frequency  interpolation  algorithms  developed  above. 


1.  Calculate  the  signal  DFT,  the  DFS. 

2.  Find  the  DFS  maximum  max(|X(A;)|)  =  |X(p)|,  and  check  if  above  threshold. 

3.  Interpolate  signal  frequency,  amplitude,  and  phase  using  interpolation  algorithms. 

4.  Subtract  estimated  spectral  leakage  of  large  signal  values  from  DFS  . 

5.  Check  if  the  maximum  subtracted  DFS  bin,  ma x.(XsLR.(k))  =  Xslr(p),  is  above 
threshold. 

For  the  above  interpolation  algorithms,  the  amplitude  of  the  sine  peak,  a,  is  estimated 
as  [20] 


a=\B\ 


(1-|A|) 


7tA 


(1  —  2|  A|  +  2A2)  ysin(7rA) J 
where  B  is  a  weighted  combination  of  the  peak  X(p)  and  the  peak’s  largest  neighbor 


(3.33) 


B  =  (l-\A\)X(p)-\A\X(p  +  a).  (3.34) 


For  real  signals,  the  signal  amplitude  is  obtained  from  the  sine  peak  amplitude  using 
Ai  =  The  signal  phase  is  estimated  as  [20] 

i fa  =  arg  (B)  —  7tA.  (3.35) 


3-4-1  SLR  Error  Analysis.  As  stated  previously,  the  effectiveness  of  the  SLR 
method  hinges  on  how  well  the  side-lobes  have  been  estimated,  which  depends  on  the 
frequency  (i.e.  centering),  amplitude,  and  phase  estimates.  Figure  3.19(a)  contains  a 
plot  of  the  mean  spectral  leakage  estimate  error  versus  bins  and  Signal  1  SNR  from  a 
1000  trial  MC  trial  using  the  GPBI  algorithm  for  frequency  estimates  for  the  following 
simulation  parameters:  [f\  =  32  +  A  bins  ,<fii  =  U[ 0,  27t]],  Signal  2  SNR=10  dB,  [f\  = 
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35  +  A  bins,  (^>2  =  C/[0, 27r]] ,  N=256.  The  A  values  are  uniformly  distributed  between 
—0.5  and  0.5  for  the  MC  simulation.  The  estimated  DFS  spectral  leakage  for  signal  1  is 
estimated  as  (the  real  system  would  use  a  look-up  table) 

Si(k)  =  DFS{ylicos(27r/i  +  </d)}  (3.36) 

The  error  is  calculated  in  dB  as 

1  M 

e%LR(k)  =  20  log10(—  £  | Sr  (A;)  -  S\(k)\)  (3.37) 

1=1 

where  e^LR(fc)  is  the  mean  error  for  the  kth  bin.  The  mean  estimated  spectral  leakage 
error  is  extremely  small  except  near  the  peak  which  is  expected  since  any  percentage  error 
is  higher  near  the  peak  since  the  spectral  leakage  is  higher  (the  peak  and  its  two  neighbors 
used  for  interpolation  are  ignored  for  detection  purposes  in  the  SLR  method).  Also  the 
mean  error  increases  as  SNR  1  increases  which  is  also  expected  because  the  magnitude  of 
the  spectral  leakage  increases. 

For  comparison  to  the  mean  error,  the  signal  1  spectral  leakage  versus  bins  and 
Signal  1  SNR  is  plotted  in  Fig.  3.19(b).  The  mean  error  is  close  to  55  dB  down  from  the 
associated  Signal  1  Spectral  leakage  in  most  places. 

Figure  3.19(c)  contains  a  plot  of  the  estimated  spectral  leakage  for  signal  1.  The 
variance  is  relatively  small,  except  near  the  signal  1  peak  which  is  expected  due  to  the 
large  swing  in  values  from  varying  A  values. 

The  above  analysis  confirms  that  the  SLR  method  compensates  for  the  signal  1 
spectral  leakage.  Thus,  the  next  step  is  to  determine  how  well  the  SLR  method  performs 
estimation/detection  of  the  signal  2  through  simulation. 

3-4-2  SLR  Simulation  Description.  With  so  many  variables  in  play,  the  best 
method  is  for  a  robust  MC  simulation  to  determine  the  Probability  of  Detection  Pci  and 
the  Probability  of  False  Alarm  Pfa  due  to  spectral  leakage  for  a  given  noise  threshold.  The 
simulation  is  setup  as  follows: 
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Figure  3.19:  Spectral  Leakage  Estimate  Error  Statitics  vs.  Signal  1  SNR  and  Bins. 

Colorbar  indicates  dB  value  of  parameter. 
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1.  Generate  M  N-length  AWGN  sequences,  where  M  is  the  number  of  MC  trials  (set  an 
arbitrary  noise  power,  in  the  plots  a2  =  0.01  is  used). 

2.  Calculate  the  noise  DFS  of  each  N-length  sequence. 

3.  Set  the  noise  threshold  by  taking  the  max  magnitude  of  each  of  the  M  N-length 
noise  sequences  for  the  positive  frequency  bins.  Then,  sorting  the  max  values  by 
magnitude  and  setting  the  threshold  at  the  10th  highest  value  (for  a  noise  only  Pfa 

<*$)■ 

4.  Specify  center  bins  for  two  signals. 

5.  Generate  M  uniformly  distributed  Ai  and  A2  values  and  add  to  the  corresponding 
center  bins. 

6.  Generate  the  M  N-length  two  sinusoids  sequences  for  the  generated  A  values  for  a 
specified  Signal  1  SNR  and  Signal  2  SNR. 

7.  Calculate  the  DFS  of  the  2  sinusoids  plus  noise  measurements. 

8.  Interpolate  the  Signal  1  frequency,  amplitude  and  phase  using  one  of  the  frequency 
interpolation  algorithms  (assuming  Signal  l’s  amplitude  is  above  threshold,  which  it 
is  for  all  simulations). 

9.  To  simulate  a  perfect  look-up  table,  calculate  the  the  signal  1  spectral  leakage  with 
estimated  frequency,  amplitude,  and  phase  as  in  (3.36). 

10.  Subtract  the  estimated  signal  1  DFS  from  the  DFS  of  the  2  sinusoids  plus  noise  and 
take  the  absolute  value  as  in  (3.20). 

11.  Perform  a  threshold  check  to  see  if  any  bin  values  exceed  the  threshold,  ignoring  the 
max  Signal  1  DFT  bin  and  its  two  closest  neighbors';  i.e.  the  ones  involved  with 
the  Signal  1  interpolation,  a.)  If  the  threshold  is  not  exceeded,  nothing  is  declared 
(Signal  2  is  not  detected),  b.)  If  the  threshold  is  exceeded  for  any  bin,  find  the 
maximum  bin  value.  If  the  maximum  bin  value  is  within  1  bin  of  the  signal  center 
bin  a  detection  is  declared,  i.e.  if  true  /2  is  40.3,  the  threshold  exceeded  at  bins  39, 

7Signals  within  2  bins  are  considered  high  resolution  in  this  thesis  and  are  not  considered  for  the  DFT 
intercept  receiver,  they  are  analyzed  in  the  IGLS  parametric  receiver  in  Chapter  V. 
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40,  and  41  would  all  be  declared  detections,  c.)  If  the  threshold  is  exceeded  by  other 

bins,  a  false  alarm  is  declared. 

3.4.3  SLR  Simulation  Results.  Figure  3.20  contains  the  Pd  and  Pfa  for  the 
MC  simulation  described  above  with  the  following  simulation  parameters:  [f\  =  32  + 
A  bins  ,  0  1  =  C/[0,27t]],  [f±  =  36  +  A  bins,  02  =  U[0,2n]],  N=256,  A  =  U[— 0.5, 0, 5].  The 
second  signal  SNR  is  fixed  at  10  dB  (for  a  sine  peak  that  is  well  above  the  detection 
threshold),  and  the  first  signal  SNR  is  increased  at  1  dB  increments.  The  threshold  for  the 
simulation  is  set  at  0.01  Pfa  for  noise  only  by  generating  the  1000  trials  of  noise,  thus  the 
false  alarms  are  for  sine  sidelobes  from  signal  1  exceeding  the  threshold.  The  second  signal 
is  detected  with  a  Pd  of  1  until  the  signal  1  SNR  exceeds  63  dB  for  the  MPP  and  PBI 
algorithms,  and  67  dB  for  the  GPBI  algorithm.  Detection  and  false  alarm  performance 
deteriorates  from  this  Signal  1  SNR  value  up.  Thus  for  a  Pd  of  1,  the  IDR  when  the  signals 
are  this  close  is  63  —  10  =  53  dB  for  the  MPP  and  PBI  algorithms,  and  67  —  10  =  57  dB 
for  the  GPBI  algorithm.  The  MPP  and  PBI  performance  is  similar  because  their  signal  1 
estimation  performance  is  similar  except  for  low  A  where  there  is  little  spectral  leakage  to 
contend  with.  The  GPBI  slightly  outperformed  the  MPP  and  PBI  algorithms  across  all  A 
values,  thus  the  GPBI  has  a  slightly  higher  IDR. 

Figure  3.21  contains  the  Pd  and  Pfa  for  the  MC  simulation  described  above  with 
signal  2  centered  at  bin  64.  The  performance  mirrors  the  performance  for  bin  36  in  3.21. 

Figure  3.22  contains  the  Pd  and  Pfa  for  the  MC  simulation  described  above  with 
signal  2  centered  at  bin  100.  The  performance  mirrors  the  performance  for  bin  36  in  3.21. 
Thus,  using  the  SLR  method,  the  IDR  has  little  dependence  on  frequency  separation,  which 
is  an  excellent  property.  This  frequency  separation  independence  can  be  can  be  explained 
by  examining  the  histogram  of  the  maximum  bin  for  the  compensated  DFT  at  an  SNR 
of  70  dB  in  Figure  3.23.  The  detections  are  clustered  around  bin  100  as  expected,  and 
the  false  alarms  are  all  clustered  at  bin  30  and  34,  which  is  where  the  spectral  leakage  is 
greatest  in  the  signal  1  DFT  (recall  that  bins  31,  32,  and  33  are  not  considered  because 
they  are  used  by  the  interpolation  algorithm).  Which  means  that  because  of  the  signal 
1  magnitude,  the  compensated  spectral  leakage  of  signal  1  is  exceeding  the  signal  2  bins 
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Table  3.1:  SLR  IDR  results  (Noise  only  Pfa  =  0.01) 


freq 

Signal  2  SNR 

Freq  Sep 

GPBI  IDR 

freq  <1.5  bins 

10  dB 

4  bins 

57  dB 

freq  <1.5  bins 

10  dB 

32  bins 

57  dB 

freq  <1-5  bins 

10  dB 

68  bins 

57  dB 

maximum  height,  most  likely  for  large  A  values.  The  Pd  and  Pfa  gracefully  degrades  as 
lower  and  lower  A  values  compensated  spectral  leakage  exceeds  the  signal  2  bins  maximum 
height.  The  false  alarms  are  clustered  around  the  same  2  bins  for  simulations  when  signal 
2  is  centered  around  bins  64  and  36  also. 

In  Figure  3.24,  the  first  signal  SNR  is  set  at  15  dB,  and  the  SNR  of  Signal  2  is  slowly 
increased  until  detected  with  a  Pd  of  one  at  -2  dB  for  simulation  parameters:  [f±  =  32  + 
A  bins  ,  (pi  =  U[ 0,27r]],  [/i  =  100  + A  bins,  =  U[ 0,27r]],  N=256,  A  =  U[— 0.5,  0,  5].  This 
can  be  considered  the  detection  threshold  limited  detection  scenario.  Notice  that  all  three 
methods  exhibit  the  same  performance  since  the  side-lobes  of  signal  are  well  compensated 
by  all  three  methods  to  below  the  detection  threshold.  This  threshold  prevents  false  alarms 
from  the  noise,  and  thus  is  not  considered  to  limit  IDR. 

3-4-4  SLR  IDR  Results.  The  SLR  method  provides  coarse  numerical  frequency 
estimates  which  is  an  EW  requirement.  Table  3.1  summarizes  the  SLR  IDR  results  for 
the  thesis  definition  of  IDR.  In  the  table,  freq  is  the  required  detection  frequency  accuracy 
(note  that  since  the  intercept  receiver  is  detecting  and  quantizing  that  this  is  not  the  Mean 
Square  Error  accuracy).  This  number  is  worst  case,  most  likely  the  worst  case  accuracy 
is  freq  <  1  bin.  Future  work  could  analyze  employing  an  interpolation  algorithm  on  the 
signal  2  peak  to  improve  frequency  estimates. 

3.5  Conclusion 

A  DFT-based  intercept  receiver  is  analyzed  without  and  then  with  noise.  With  no 
noise,  it  is  shown  that  the  DFT-based  intercept  receiver  IDR  is  limited  by  spectral  leakage 
that  is  exclusively  a  result  of  finite  measurement  time.  To  handle  spectral  leakage,  the  novel 
SLR  method  is  analyzed  for  its  IDR  performance  in  white  noise  using  a  detection  threshold 
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Figure  3.23:  Histogram  of  SLR  compensated  DFT  detections  (Noise  alone  Pfa  =  0.01, 
Signal  1  bin  32,  Signal  2  bin  100). 


Figure  3.24:  SLR  method  Pd  for  low  SNR  1  (Noise  alone  Pfa  =  0.01,  Signal  1  bin  32). 
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scheme.  The  method’s  IDR  performance  is  shown  to  have  little  frequency  separation 
dependence  and  provides  high  IDR  with  the  required  numerical  frequency  estimates. 

Because  of  the  spectral  leakage  biasing  the  estimates,  it  is  very  difficult,  if  not  impos¬ 
sible,  to  achieve  optimal  (ML)  estimates  using  only  the  DFT.  However,  if  the  number  of 
sinusoids  present  in  the  measurement  is  known  already,  optimal  estimates  can  be  obtained. 
This  optimal  estimate  analysis  is  the  topic  of  the  next  two  chapters. 
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IV.  Cramer- Rao  Bound  for  Instantaneous  Dynamic  Range 


!  For  Chapters  IV  and  V  in  the  thesis,  the  number  of  signals  present  in  the 

measurement  is  assumed  known. 

4.1  Introduction 

Before  the  discussion  of  the  IGLS  algorithm-based  parametric  receiver  in  Chapter 
V,  a  natural  question  is  what  is  the  best  IDR  that  can  be  achieved  for  a  desired  RMS 
frequency  estimation  accuracy.  For  unbiased  estimators,  the  Cramer-Rao  Bound  (CRB)  is 
used  to  provide  a  lower  bound  on  the  MSE  of  the  estimates  [11-13].  The  CRB  is  a  lower 
bound  for  the  variance  of  any  unbiased  ML  Estimator.  Thus,  the  IDR-CRB  is  derived 
below.  Section  4.2  derives  the  multiple  frequency  CRB  originally  derived  by  Rife  [6].  The 
multiple  frequency  CRB  is  modified  by  the  author  in  Section  4.3  to  arrive  at  the  IDR-CRB 
for  complex  and  real  signals. 

4-2  Derivation  of  the  CRB  for  Multiple  Sinusoids  in  AWGN 

In  this  section,  the  CRB  for  multiple  sinusoids  in  white  noise,  originally  derived 
by  Rife  in  [6],  is  derived.  First,  in  Section  4.2.1,  the  multiple  sinusoid  CRB  is  derived 
for  complex  signals.  Next,  in  Section  4.2.2,  the  complex  signals’  CRB  [6]  is  extended  to 
develop  the  real  signals’  CRB. 

4.2.1  Complex  Signal  CRB.  Consider  the  following  complex  sinusoidal  signal 
vector,  sc,  with  unknown  parameters 

p  p 

sc(n)  =  ^  ^exp{  j{ojjn  +  (f)j)}  =  n  =  0  . . .  N  -  1  (4.1) 

i= 1  2—1 

Let  the  unknown  parameter  vector  6  be  defined  as 

0  =  [ui  Alfa...  up  Ap  fa\ .  (4.2) 
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The  signal  sc  is  combined  in  a  transmission  channel  with  noise  to  form  the  measurement 
vector 


xc  =  sc  T  wc,  (4.3) 

where  wc  is  complex  noise  distributed  as 


wc  =  Af(0,  R). 


(4.4) 


Since  xc  is  a  linear  combination  of  a  Gaussian  distributed  random  variable,  xc  is  itself  a 
Gaussian  distributed  random  variable 


xc=A7(sc,R).  (4.5) 

Thus,  the  Probability  Density  Function  (PDF)  of  xc  (with  the  complex  noise  assumption) 
is  [13] 

fo(xc)  =  (vrcr2)^Ar|R|_1exp  {— (xc  -  sc)//R_1(xc  -  sc)}  .  (4.6) 

Using  the  AWGN  assumption,  Equation  (4.6)  simplifies  to 

/e(xc)  =  (vrcr2)^Arexp  |--^(xc  -  sc)^(xc  -  sc)|  .  (4.7) 

For  the  Cramer-Rao  development,  the  natural  log  of  (4.7)  is  helpful, 

In  (/0(xc))  =  Le{xc)  =  -Mn(vrcr2)  -  (xc  -  sc)H  (xc  -  sc),  (4.8) 

where  Lg(xc)  is  called  the  log  likelihood  function  (referred  to  as  Likelihood  because  the 
value  of  the  unknown  parameter  vector  6,  which  characterizes  the  transmitted  signal  sc, 
is  estimated  by  determining  the  value  that  made  the  known  measurement  vector  xc  most 
likely  to  occur  [10]). 

The  derivative  of  (4.8)  is 


s(0,xc) 


dL(0,xc) 

86 


(4.9) 
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where  s(0.  xc)  is  called  the  score  function.  The  values  of  6  where  the  score  function  vanishes 
is  the  ML  estimate  of  6  -  see  Section  2.4.1.  The  score  function  covariance  is  referred  to  as 
the  Fisher  information  matrix,  J(0) 

J(0)  =  £  {s(6,xc)s(0,xc)H}  .  (4.10) 


A  per-element  formula  for  the  Fisher  information  matrix  is  [11] 


JiPuOj) 


f  O2L(0,xc)  } 

I  09,09 j  J ' 


(4.11) 


The  Fisher  Information  contains  information  on  the  maximum  rate  of  change  near  the  peak 
of  the  pdf  which  corresponds  to  the  ML  estimate.  The  inverse  of  the  maximum  rate  of 
change  yields  the  minimum  variance  the  unbiased  estimate  can  attain,  thus  the  inverse  of 
the  Fisher  information  matrix  contains  information  on  the  minimum  value  the  covariance 
can  attain  [10].  If  the  estimate  is  unbiased,  the  inverse  is  also  the  minimum  MSE  a  ML 
estimate  can  attain  [10]  -  see  Section  2.3. 

A  more  specific  formula  for  the  per  element  Fisher  information  matrix  of  xc  is  derived 
below.  Inserting  (4.8)  into  (4.11)  yields 


Jx[Pi,0j) 


(4.12) 


Using  the  chain  rule  from  Calculus,  (4.12)  becomes 


Me,M^£,^iir^)iKc-Sc)  +  (Xc-Scf(,d^  .  ■  >0-  ■ 


1  „  02 s? 


c72yyOOiOOj 


'  06,09 j 


)  +  (^)(^)  +  (ii2r)(W:))  • 


06i  09 


09 j  09i 


(4.13) 


Taking  the  expectation  of  the  terms  of  (4.13)  and  noting  that  following  expectation  value, 
£  {xc  —  sc}  =  0,  the  formula  for  the  Fisher  information  matrix  simplifies  to 


1  „ ,  <9s^  ,  ,  Os, 


Ost1 ,  ,  0sr 


00 j  09i ' 


(4.14) 


Let 


0sc 

09 i 


a  +  jb 
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(4.15) 


where  a  is  the  real  part  and  b  is  the  imaginary  part  of  the  complex  vector  quantity.  Also, 


let 

Wi=C  +  jd  (4.16) 

where  c  is  the  real  part  and  d  is  the  imaginary  part  of  the  vector  quantity.  Then  (4.14) 
becomes 

Jx(0i,  Qj)  =  ^ ((a  +  jb)H(c  +  jd)  +  (c  +  jd)H(a  +  jb)).  (4.17) 

Equation  (4.17)  simplifies  to 

Jx{0i,  Oj)  =  ^(aTc  +  dTb).  (4.18) 


Inserting  the  derivatives  from  (4.15)  and  (4.16)  into  (4.18),  the  desired  formula  for  the 
Fisher  information  matrix  for  complex  sinusoids  in  complex  white  Gaussian  noise  is  ob¬ 
tained 


Jx(0i,  Qj)  =  ~^Re 


a z 


ds^  dsc 
ddj  ddj 


=  -^ReiSj} 


u 


(4.19) 


where 

=  gsf  dsc 
3  86i  89 j  ' 


(4.20) 


The  above  formula  is  used  to  calculate  the  Fisher  information  matrix  elements  for 
the  following  9  vector  parameters: 


sci  SC2 


0(1)  =  U\  0(4)  =  L02 
0(2)  =  Ai  0(5)  =  A2 
0(3)  =  0i  0(6)  =  02 
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The  hermitian  matrix  S?  contains  the  values 


4  2dhf  Oh, 

■^*■1  duJi  duj\ 

a  dh? , 

a2&^\  <9hi 

1  du>i  dip i 

A  .  A®*1?  Oh-2 
^lji2  dun  duJ2 

a  dh? , 

auh2 

9ha 

1  2  9^  902 

(1,2)* 

hfhr 

A  i  Oh!? 

hfh2 

x  i  dh? 

(1,3)* 

(2,3)* 

A  2  Oh?  Oh! 

1  dipi  dipi 

A  A  Oh?  9h2 
dip i  duj2 

a  Oh? , 

^iwr*12 

2  1  dipi  dtp2 

(1,4)* 

(2,4)* 

(3,4)* 

A  2  Oh?  0h2 

2  du2  Olo 2 

a  dh? , 
A2~a^h2 

A  2  9hf  9h2 

2  du}2  dip2 

(1,5)* 

(2,5)* 

(3,5)* 

(4,5)* 

hfh2 

M  hfij 

(1,6)* 

(2,6)* 

(3,6)* 

(4,6)* 

(5,6)* 

/)2  9h|f  9h2 
^2  c?<^>2  902 

l.  y  4  -r  4  -i 

(4.21) 

The  structure  of  S?-  allows  the  Fisher  information  matrix  (4.19)  to  be  factored  as  [6] 


J  =  -^CMC, 

C 7Z 

where  C  is  the  diagonal  matrix  of  amplitudes 


_ 

_ 

Ai 

0 

0 

Cl 

0 

c  = 

,  Cj  = 

0 

1 

0 

0 

c2_ 

L 

0 

0 

Ai 

(4.22) 


(4.23) 


The  M  matrix  from  (4.22)  incorporates  the  h  terms  of  (4.19)  and  has  the  following  block 


form 


M  = 


Mu 

M2i 


M12 

M22 


(4.24) 


Mi;j  is  a  3  by  3  matrix  that  is  derived  in  the  following  manner.  Note  the  following 


dh^did^^  =  ^nexP{i(w*n  +  fa)}  (4.25) 

dh^  —  =  iexP  {j{uin  +  fa)}  (4.26) 

ocpi 

Using  (4.25),  the  elements  of  are 

ou  H  ou  .  N~1 

Mij(  1,1)  =  i?e{— - —  — — }  =  V  ?r2cos((wj  -Uj)n  +  (<f>i  -  <f>j))  (4.27) 

dui  ocoj  n=0 
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-  <9h 


H 


N- 1 


My (1,2)  =  Re{-±-h.j}  =  ^2  -nsin((ui  - uj)n+  (fa  ~  fa))} 


duJi 

-  dh/  fill 


n= 0 
JV-1 


My  (1,3)  =  Re{— —  — ^  ncos((wj  -  wj)n  +  (fa  -  ^)) 


<9(/> 

r  9h 


»i=0 

JV-1 


My  (2, 1)  =  i?e{hf ^  nsin((wi  -  ujfan  +  (fa  -  <£,■)) 

•?  n=0 

JV-1 

My (2, 2)  =  i?e{hfhj}  =  ^  cos ((cu*  -  Uj)n  +  (fa  -  fa)) 


n= 0 
JV-1 


My(2,  3)  =  i?e{hf  — ^}  =  ^  sin((o;i  -  Uj)n  +  (fa  -  fa)) 

Vi  n= 0 

ou-fj  fiU  , 

My  (3, 1)  =  i?e{— =  V  ncos((aii  -  Lofan  +  (4  -  </»_,-)) 
dw,  n=0 

<9hff  iV_1 

My  (3, 2)  =  i?e{— -Mij}  =  X]  _sin((wi  ~  Wj)n  +  (0i  -  <£,-)) 


dfa 

,  dh.1/  dh 


n= 0 
JV-1 


My (3, 3)  =  ^  cos(((jj  -wi)n+  (fa  -  ^)) 


n=0 


Let  Ay  =  (oij  —  ojj)n  +  (fa  —  4>j),  then  the  matrix  form  of  My  is 


My  — 


Etc1  ™2cos(Ay)  Etc1  -nsin(Ay)  En=o  ^cos(Ay) 
tJV-1 


EnEo1  nsin(Ay)  E^o  oos(Ay)  En=o  sin(Ay) 
Y)n=o  ncos(Ay)  Et~o  —sin (Ay)  EnE?  cos(Ay) 


(4.28) 

(4.29) 

(4.30) 

(4.31) 

(4.32) 

(4.33) 

(4.34) 

(4.35) 

(4.36) 


The  inverse  of  the  Fisher  Information  Matrix  is  the  minimum  estimation  error  co- 


variance  matrix 

J1  =  — C-1M_1C-1. 
2 


(4.37) 
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Each  element  of  the  diagonal  of  (4.37)  is  the  CRB  for  the  corresponding  parameter.  It  can 
be  shown  that  a  formula  for  the  inverse  of  the  blocked  matrix  M  is  [13] 


M_1  = 


(Mn  - 


-(Mu  -  Mi2M221M2i)-1M12M221 


— (M22  —  M2iM111Mi2)  1M2iM111  (M22  -  M2iM111M 


12) 


\-i 


(4.38) 

Due  to  the  symmetry  of  M  for  the  complex  sinusoid  case,  M_1(l,l)  =  M_1(4, 4)  [6]. 
Thus,  the  CRB  for  minimum  u>i  estimation  variance  is 


varjchi}  > 


v2P(l,l) 

2  A2 


(4.39) 


where  D  is  defined  as 


D  —  (Mu  —  Mi2M221M2i) 


-i 


(4.40) 


The  corresponding  CRB  for  frequency  estimation  is 

°2D(  1,1) 


var {fi}  > 


2(2vr)2Af  ■ 


(4.41) 


To  ensure  the  computer  code  generating  the  multiple  sinusoid  CRB  accuracy,  com¬ 
parisons  to  the  results  in  [13]  and  [6]  are  contained  in  Figure  4.1.  Figure  4.1(a)  contains  the 
CRB  versus  SNR  for  two  complex  sinusoids  in  complex  noise  with  the  following  parame¬ 
ters:  N=25,  [^i  =  l,/i  =  0.5 ,</>!  =  0],  [A2  =  l,/2  =  0.52,  ^>2  =  f];  SNR  =  -10  logger2). 
The  results  match  those  in  [13].  Figure  4.1(a)  contains  the  CRB  versus  frequency  for  one 
real  sinusoid  with  worst  phase  difference  between  negative  and  positive  frequency  images 
(which  consists  of  two  complex  sinusoids))  in  real  noise  with  the  following  parameters: 
[A\  =  1,  0i  =  0,  0ij  =  Nirf-2  —  ir  —  Nirfi]-,  SNR,  =  —  101ogi0(2er2)  =  20.  Results  match  [6]. 
Thus,  the  code  generating  the  multiple  sinusoid  CRB  is  validated. 


4-2.2  Real  Signal  CRB.  For  real  signals  in  real  white  noise,  the  development 
from  above  is  modified.  The  signal  is  now 

p  p 

s(n)  =  T.  Ajcosjujt  +  (j)j)  =  Ajg(uj,  (j)j ).  n  =  0 . . .  IV  —  1  (4.42) 

i—  1  i=  1 
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Equation  (4.22)  becomes  [6] 


J  =  -icQC 

ga 


(4.43) 


where  Q  has  the  same  block  structure  as  (4.24).  Each  block  matrix  of  elements  of  Q  is 
defined  as  [6] 

Q ij  =  2  4"  ^ji  4*i  4“  4>j)^j\i  (4.44) 

where  was  defined  in  (4.36),  and  Bj  is  a  3  by  3  diagonal  matrix  defined  as 


B  j  =  diag(l,  —1, 1). 


(4.45) 


Because  the  tones  are  real  instead  of  complex,  J_1(4,4)  is  not  in  general  equal  to  J-1(  1, 1) 
(usually  only  a  slight  difference).  Thus,  a  formula  for  each  frequency  is  used  for  the  bound 
on  the  corresponding  frequency  estimate  of  real  signals.  The  Cramer-Rao  bound  for  the 
first  frequency  estimate  is 


var{/i}  > 


a^Ej  1,1) 

(2vr)M? 


(4.46) 


where  E  is  defined  by  substituting  Q)?  for  My  in  (4.40).  The  Cramer-Rao  bound  for  the 
second  frequency  estimate  is 


var{/2}  > 


v2F(l,l) 

(2vr)2Mi 


(4.47) 


where  the  matrix  F  is  defined  as 


F  =  (Q22  —  Q2iQn1Qi2)  1- 


(4.48) 


4-3  Cramer-Rao  Bound  for  Instantaneous  Dynamic  Range 

In  [5],  the  IDR-CRB  for  complex  signals  is  derived.  The  method  of  [5]  uses  an 
iterative  method  in  terms  of  delta  values  for  the  calculation  of  the  IDR-CRB,  which  is 
computationally  intensive  and  difficult  to  implement.  Because  the  Fisher  information 
matrix  can  be  factored  as  shown  in  (4.22),  a  simpler  and  more  reliable  method  is  introduced 
here  to  calculate  the  complex  signal  IDR-CRB  and  then  extended  to  the  real  signal  IDR- 
CRB. 
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Recall  the  definition  of  IDR  employed  in  this  thesis 


•  The  thesis  IDR  definition  -  IDR  is  defined  as  the  maximum  signal  difference  for 
a  given  frequency  estimation  accuracy,  a  given  frequency  separation  and  a  given 
SNR  [5]. 

From  this  definition,  the  CRB  for  instantaneous  dynamic  range  for  two  complex  exponen¬ 
tials  is  derived  by  modifying  the  multiple  sinusoid  CRB  using  the  following  method. 


1.  Specify  the  noise  power,  a 2,  and  the  desired  SNR  (in  terms  of  A\),  RMS  frequency 
estimation  accuracy1,  and  frequency  separation  A/. 

2.  The  A]  amplitude  is 

Ai  =  \J  <7210(5jvr/10).  (4.49) 


3.  Use  the  following  modified  form  of  (4.41)  to  solve  for  the  bound  amplitude  of  A 2 


A-2b  > 


°2D(  1,1) 

2(27 TYPaa 


(4.50) 


where  facc  is  the  desired  RMS  frequency  estimation  accuracy  and  A2b  is  the  bound 
amplitude.  If  A2b  is  greater  than  A\,  then  the  desired  frequency  accuracy  is  not 
achievable  for  the  given  parameters. 

4.  If  achievable,  the  Cramer  Rao  bound  for  the  instantaneous  dynamic  range  is  then 
defined  as 

IDR(dB)<201oglo(-^-)  .  (4.51) 


Figure  4.2  is  a  plot  of  the  complex  signal  IDR-CRB  vs.  frequency  difference  for 
the  parameters:  [Ai  =  l,/i  =  0.2,  fa  =  0],  [<j)2  =  Nirf2  -7 r  -  ADr/i];  facc  = 

SNR  =  — 101og10(cr2)  =  20,  N  =  100.  The  phase  value  f>2  used  to  generate  Fig.  4.2 


<t>2  =  Nnf2  -ir  -  Nirfi, 


(4.52) 


1Recall  from  Chapter  2.3  that  the  RMS  squared  is  the  MSE,  which  is  equal  to  the  estimate  variance  for 
an  unbiased  estimator. 
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30 


Figure  4.2:  Complex  IDR-CRB  versus  A /  . 


is  the  worst  phase  difference  for  frequency  estimation  [5].  The  point  where  the  plot  flat¬ 
tens  out  where  the  phase  difference  and  frequency  difference  interaction  terms  M  become 
negligible  and  the  single  signal  required  amplitude  is  attained  for  the  desired  estimation 
accuracy. 

The  method  of  calculating  the  real  signal  IDR-CRB  is  similar  to  the  method  for 
complex  signals,  with  the  following  exceptions: 

•  The  Aj  is  now  calculated  as 


Ai  =  v/2cr210('SAri?/1°). 


(4.53) 


The  following  equation  must  be  used  in  place  of  (4.50)  to  calculate  the  bound  am¬ 
plitude 

(4.54) 


(2n  )2/', 


2  f  2 
acc 
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Figure  4.3:  Real  Signal  IDR-CRB  Algorithm  Flowchart. 


A  slight  error  term  is  introduced  using  the  above  equation  when  the  amplitude  is  near 
the  instantaneous  dynamic  range  resolvable  threshold  because  of  the  slight  difference 
in  the  values  of  J_1(l,l)  and  J-1(4,4)  for  real  tones.  The  algorithm  is  modified 
without  iterating  to  handle  this  slight  difference  in  the  following  way. 

•  Compare  the  estimation  accuracy  of  (4.46)  to  the  desired  estimation  accuracy  if  the 
value  of  A2b  is  close  to  the  value  of  A\ .  If  the  estimation  accuracy  is  higher,  the 
threshold  has  not  been  met  even  though  A 2  amplitude  is  less  than  one,  and  the 
desired  estimation  accuracy  is  not  achievable. 

•  The  IDR-CRB  values  for  real  signals  are  frequency  dependent,  especially  at  frequen¬ 
cies  near  zero  (i.e.  if  f\  =  O.OOlj. 

Figure  4.3  is  a  flow  chart  representation  of  the  method  to  generate  the  IDR-CRB. 

!  Throughout  the  rest  of  the  thesis,  IDR-CRB  denotes  the  real  signal  IDR-CRB 
unless  otherwise  stated. 

Figure  4.4  is  a  plot  of  the  real  signal  IDR-CRB  for  different  values  of  N  generated 
using  the  method  of  Figure  4.3.  In  the  plots,  the  IDR-CRB  is  calculated  at  A /  =  0.001 
intervals  for  specified  N  samples  of  two  sinusoids  in  noise  with  worst  phase  difference  with 
parameters:  [A\  =  l,fi  =  0.2, (j>i  =  0],  [<^>2  =  ADr/2  —  7r  —  ADr/i];  SNR  =  — 10 log10(2(j2)  = 
20  dB.  For  comparison  purposes,  the  frequency  estimation  accuracy  is  also  scaled  as  a 
function  of  N,  where  N  is  the  number  of  measurement  points.  As  expected,  the  bound  for 
the  higher  N  exhibits  a  sharper  rise  time  than  the  lower  N.  In  the  Chapter  V,  an  IGLS 
algorithm-based  parametric  receiver  IDR-CRB  comparison  validates  the  IDR-CRB  results. 
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4-4  Conclusion 

A  Cramer-Rao  bound  is  developed  for  an  unbiased  estimator’s  IDR  by  modifying 
Rife’s  result  in  [6]  for  both  real  and  complex  sinuoidal  signals  in  AWGN.  To  achieve  the 
IDR-CRB,  the  estimator  must  be  an  unbiased  ML  estimator  of  frequency;  i.e.  the  estimator 
must  be  unbiased  and  efficient  [11].  The  next  chapter  compares  an  IGLS-algorithm  based 
parametric  receiver  to  the  IDR-CRB  derived  in  this  chapter. 
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V.  IGLS  Algorithm-Based  Parametric  Receiver 

5. 1  Introduction 

When  the  number  of  signals  is  known,  the  EW  receiver  frequency  estimates  are  im¬ 
proved  dramatically  by  basing  the  receiver  on  a  parametric  frequency  estimation  algorithm, 
i.e.  a  parametric  EW  receiver  1.  The  parametric  algorithms  achieve  these  improved  esti¬ 
mates  because  prior  knowledge  of  the  signal  form  is  exploited.  To  bound  IDR  performance 
and  compare  to  the  IDR-CRB  derived  in  Chapter  IV,  a  parametric  frequency  estimation 
algorithm  that  achieves  ML  results  is  desired.  Thus,  a  frequency  estimation  algorithm 
based  on  Linear  Prediction  (LP)  called  Iterative  Generalized  Least  Squares  is  introduced 
for  application  in  the  parametric  receiver  and  shown  to  yield  Maximum  Likelihood  fre¬ 
quency  estimates  in  Section  5.2.  The  IGLS  algorithm  is  compared  to  the  IDR-CRB  in 
Section  5.3.  Finally,  in  Section  5.3.1,  experimental  results  from  the  IDR-CRB  compari¬ 
son  result  in  defining  IDR  differently  for  a  Parametric  based  receiver  when  the  frequency 
estimate  requirements  are  loose. 

5.2  IGLS  Development 

In  this  section,  the  IGLS  algorithm,  related  to  the  IQML  algorithm  discussed  in 
[14],  is  fully  developed  and  shown  to  yield  ML  frequency  estimates.  This  is  the  author’s 
development  of  the  IGLS  algorithm  originally  derived  by  Dr.  Pachter  and  researched 
by  Ingham  in  [8]  and  Zahirniak  in  [7].  Section  5.2.1  provides  the  necessary  background 
on  LP  theory  including  Prony’s  method  and  the  Extended  Prony  Method.  Section  5.2.2 
derives  the  IGLS  algorithm.  In  Section  5.2.3,  frequency  estimate  confidence  intervals  are 
developed.  Simulations  then  verify  IGLS  ML  performance  in  Section  5.2.4. 

5.2.1  Linear  Prediction  Theory.  Consider  the  real  sinusoidal  signal 

p 

s(n)  =  Y,  Aj  cos {cjjn  +  4>j ) .  (5.1) 

i= 1 

1  Determining  the  number  of  signals  for  a  parametric  receiver  is  an  area  of  research  in  its  own  right  and 
is  not  discussed  here.  See  Reference  [13]  for  more  details. 
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The  signal’s  samples  satisfy  the  difference  equation  [21,22] 

2  p 

s(n)  =  J>(m)a(ra-m),  (5-2) 

771=  1 

where  the  coefficients  of  the  2p  length  vector  a  are  called  the  LP  coefficients.  The  proof 
for  this  relationship  is  quite  involved,  but  very  informative  (this  proof  follows  [21]  closely 
with  an  example  at  the  end  to  clarify). 

Consider  now  the  following  factored  polynomial,  3>(z),  with  the  real  signal  frequencies 
of  (5.2)  as  roots  (for  a  complex  signal  remove  the  negative  exponential), 

p 

$(z)  =  JJ(z  -  e>Ui)(z  -  (5.3) 

i— 1 

Expanding  the  above  equation  yields  the  following  polynomial 

2  P 

$(z)  =  Y  a{m)z2p-m.  (5.4) 

777=0 

where  a(0)  is  constrained  to  be  1.  Form  a  linear  difference  equation  by  multiplying  a{m ) 
by  s(n  —  m)  and  summing  over  m  to  yield 

2  p  2  p  2  p 

Y  a(m)s(n  -  m)  =  Y  a(m)  y  (e^n“m) e?*  +  (5.5) 

777=0  777=0  7=1 

Switching  the  order  of  summations  and  making  the  following  substitution  n  —  m  =  n  — 
2p  +  2p  —  m  (realizing  that  n  >  2 p) ,  yields 


2  p 


2  p  2  p 


777=0 


Y  a{m)s{n  -m)  =  Y  -^eju)^n~2p)  ej^  Y  «( 


i= 1 


e  a'm>e 
m= 0 
2  V 


ju>i(2p-m) 


+E 


A 


2  p 


e  -r,e  jy-i  ^2  a(m)e 


-jui(2p-m) 


(5.6) 


t=i 


771—0 
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Note  that  (This  is  beautiful!!) 


2  v 

a(m)(ejuJi{2p~m))  =  0  (5.7) 

m= 0 
2  V 

J2  a{m)e-juJ^2p-m))  =  0.  (5.8) 

m= 0 

because  the  exponentials  are  roots  of  the  polynomial  from  (5.4)  as  shown  in  (5.3).  Thus, 
the  LP  equation 

2  v 

s(n)  =  —  ^2  a(m)s(n  —  m),  (5.9) 

m= 1 

is  arrived  at  (recall  a(0)  is  constrained  to  be  one).  In  matrix  format,  the  exactly  determined 
system  to  solve  for  the  polynomial  values  is 


s(2p  —  1)  s(2p  —  2)  ...  s(0) 

a(l) 

s(2p) 

s(2p)  s(2p  —  1)  ...  s(l) 

a(2) 

=  - 

s(2p  +  1) 

s(4p  —  2)  s(4p  —  3)  ...  s(2p  —  1) 

1 

© 

i _ 

s(4p  —  1) 

where  the  matrix  of  signal  values  is  2 p  by  2 p.  After  solving  for  the  2 p  LP  coefficients,  the 
frequencies  are  obtained  by  rooting  the  polynomial  (5.4)  formed  by  the  LP  coefficients. 
Note  that  the  nonlinearity  of  estimating  the  frequencies  has  been  compressed  into  the 
rooting  of  the  polynomial  comprised  of  the  LP  coefficients  a  [21].  The  above  math  can 
be  interpreted  in  the  following  way:  Equation  (5.9)  is  a  linear  difference  equation  with 
associated  characteristic  equation  (5.4)  that  has  the  homogeneous  solution  given  by  (5.3) 
[21].  The  above  LP-based  method  of  determining  frequencies  via  rooting  the  LP  coefficient 
characteristic  polynomial  is  called  Prony’s  Method  [21]2. 

Formulas  for  the  relationships  between  the  LP  coefficients  a  and  the  frequencies  f 
are  derived  below  for  the  case  of  two  real  sinusoids.  For  two  real  sinusoids,  Equation  (5.3) 

2Originally  developed  by  Gaspard  Riche,  Baron  de  Prony  in  1795  in  his  study  of  the  expansion  of  various 
gases  [21]. 
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Table  5.1:  LP  Coefficients  to  Frequency  Relationship. 


LP  Coeff 

Formula 

a(0) 

1 

a(l) 

— 2cos(aq)  —  2cos(a;2) 

a(2) 

2  +  4cos(wi)  cos(w2) 

a(3) 

— 2cos(a;i)  —  2cos(a;2) 

a(4) 

1 

becomes 


$(z)  =  (z-  ejuJl)(z  -  -  ejuJ*)(z  -  e~j^).  (5.11) 


Combining  the  two  like  frequency  terms  yields 


$(z)  =  ( z 2  —  2  cos(oq)2  +  1  )(z2  —  2  cos^)^  +  !)•  (5-12) 


Multiplying  out  (5.12)  yields  the  polynomial 


$ (z)  =  z4  —  ( 2  cos (uq ) + 2  cos (uj2 ) ) z3  +  (2 + 4 cos (oq )  cos (u>2 ))z2  —  (2  cos (wi ) + 2  cos ( L02 ) ) z + 1 . 

(5.13) 

Table  5.1  relates  the  coefficients  of  the  polynomial  to  the  frequencies.  Note  that  the  Table 
5.1  polynomial  terms  are  symmetric;  it  can  be  shown  in  general  that  the  polynomial  terms 
for  real  sinusoids  are  symmetric  [22]. 


For  two  sinusoids  a  closed  form  solution  for  the  cosines  in  Table  5.1  in  terms  of  a(l) 
and  a(2)  using  the  formulas  of  Table  5.1  can  be  obtained.  Solving  the  formula  for  a(2)  in 
terms  of  cos(u;2)  yields 


cos  (0J2)  = 


o(2)  -  2 


4cos(aq) 

Insert  the  value  for  cos(a;2)  into  the  formula  for  a(l)  to  obtain 


(5.14) 


—a(  1)  =  2cos(uq)  + 


9/  \  ttf  1 ")  ,  . 

0  =  cos  (uq)  H - — -  cos(o;i)  + 


o(2)  -  2 
2  cos(tJi) 

o(2)  -  2 


(5.15) 

(5.16) 
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The  above  is  a  quadratic  equation  with  roots  at  cos(cji)  and  at  cos(u;2)  (switch  around 
(5.14)  to  see  the  second  root).  Thus,  a  quadratic  must  be  solved,  and  the  following  rela¬ 
tionship  is  arrived  at 


cos(iui),  cos(cu2) 


— a(l)  ±  v'a(l)2  -  4a(2)  +  8 
4 


(5.17) 


When  noise  is  not  present,  the  LP  coefficients  are  estimated  perfectly  from  the  data  by 
solving  (5.10)  for  the  vector  a,  and  thereby  the  frequencies  can  be  calculated  error-free  from 
(5.17)  regardless  of  frequency  spacing,  i.e.  no  Fourier  Resolution  limit.  Conceptually,  the 
difference  between  the  Periodogram  method  and  Prony’s  method  is  that  the  Periodogram 
evaluates  certain  frequencies,  while  Prony’s  method  estimates  the  frequencies  exactly  from 
the  data  [21]. 

When  measurement  noise  is  added  to  the  system,  error  is  introduced  into  the  above 
linear  prediction  relationship  and  Prony’s  method  performs  poorly.  Therefore,  the  above 
equations  are  modified  to  handle  the  presence  of  noise. 


5.2. 1.1  Extended  Prony  Method.  The  measured  signal  with  noise  is 


x(n)  =  s(n)  +  w(n ),  n  =  0, . .  N  —  1 


(5.18) 


where  s(n)  is  defined  in  (5.1)  and  w(n)  is  AWGN.  Normally,  more  samples  are  present  than 
required  for  the  exact  solution  of  the  LP  coefficients  using  Prony’s  Method.  In  the  presence 
of  noise,  these  additional  samples  can  be  exploited  to  obtain  a  least  squares  solution  that 
washes  out  the  error  introduced  by  the  measurement  noise.  For  the  following  development, 
reshape  the  measurement  vector  x  of  (5.18)  into  a  M  —  2p  by  2p  +  1  observation  matrix  X 


x(M  —  1) 

x(M  —  2)  . . 

x(M  —  2  p) 

x= 

x(M  —  2) 

x(M  —  3)  . . 

.  x(M  —  2p  —  1) 

x(2p) 

x(2 p  —  1)  . . 

x(0) 

(5.19) 
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The  following  2p  +  1  by  p  +  1  constraint  matrix  B  is  introduced  to  exploit  the  real  signal 
LP  coefficients’  symmetry  and  reduce  the  parameter  vector  size 


B  = 


0P 
1  , 


IBp  Op 


(5.20) 


where  IBp  is  the  p  by  p  ‘backwards’  identity  matrix  defined  by  IBP  =  S(P  +  1  —  i  —  j), 
Ip  is  a  p  by  p  identity  matrix,  and  0P  is  a  p  by  1  vector  of  zeroes.  Using  B,  the  linear 
prediction  vector  a  =  [a(0)...a(2p)]T  can  be  formed  in  the  following  way 


a  =  Be*.  (5.21) 

where  the  reduced  parameter  vector  a  is  defined  as  a  =  [1  a(l) . . .  a(p)]7  .  The  observation 
matrix  X  is  multiplied  by  the  constraint  matrix  B  yielding 


Xc  =  XB  =  [x0|X0], 


(5.22) 


where  the  matrix  Xc  is  the  constrained  data  matrix,  xQ  contains  the  first  column  of  Xc, 
and  X0  contains  the  remaining  columns  of  Xc.  Using  (5.22),  the  linear  prediction  model 
is 

g0  —  X0  T  X0a0.  (5.23) 

where  the  vector  aD  is  the  vector  of  LP  coefficients,  ac  =  [a(l)...a(p)]T ,  and  the  vector  eG 
is  the  prediction  error  vector  due  to  the  noise. 

The  error  power,  ||e0|||,  is  equal  to 

|  |e0|  ||  =  (x0  +  X0a0)T(x0  +  X0a0).  (5.24) 


A  good  estimate  of  the  unknown  LP  coefficients  are  the  LP  coefficients  that  minimize  the 
expression  in  (5.24),  i.e.  the  Least  Squares  (LS)  estimate 

a0  =  —  (X^X0)-1Xj  xG.  (5.25) 
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Equation  (5.25)  is  referred  to  as  the  Extended  Prony  Method  (EPM)  [21].  Although  an 
improvement  over  the  exactly  determined  method  -e.g.  (5.10)-,  the  EPM  does  not  provide 
efficient  LP  estimates,  and  performs  poorly  in  low  SNR  [21]  as  shown  in  the  simulations 
in  Section  5.2.4. 

5.2.2  Iterative  Generalized  Least  Squares.  To  improve  upon  (5.25),  it  is  first 
necessary  to  understand  why  efficiency  is  not  achieved.  To  this  end,  the  M  by  M  -  2p 
Toeplitz  matrix  Ar  defined  by  the  LP  coefficients  is  introduced 

A1  =  Toeplitz(l,  a(l) . . .  a(2p),  0  . . .  0).  (5.26) 

Using  (5.26),  an  equivalent  representation  of  the  LP  equation  error  e0  of  (5.23)  is 

eG  =  A 1  x 

=  Ats  +  ATn  (5.27) 

=  ATn. 

Thus,  the  noise  vector  n  of  the  measurement  vector  x  is  subjected  to  a  moving  average 
process  that  yields  eQ.  The  vector  eG  is  a  colored,  normally  distributed,  zero  mean  random 
vector  with  covariance  matrix 

Ceo  =  ATCnA,  (5.28) 

where  Ceo  is  the  error  vector  covariance,  e0,  and  Cn  =  cr2I  is  the  noise  vector  covari¬ 
ance.  The  EPM  least  squares  estimate  in  (5.25)  weights  each  term’s  contribution  equally. 
However,  since  a  moving  average  process  has  been  applied  to  the  noise,  coloring  the  error 
covariance  matrix,  this  assumption  is  invalid  [7];  unfortunately,  this  is  often  overlooked  in 
system  identification  work  [23,24]. 

To  account  for  the  colored  noise,  the  extended  Prony  method  is  modified  in  the  follow¬ 
ing  way.  Perform  a  Cholesky  decomposition  of  the  positive  semi-definite  error  covariance 
matrix  inverse  to  obtain 

C-1  =  GGt,  (5.29) 
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where  G  is  the  Cholesky  decomposition  of  Ce  G  =  Ceo2 .  Equation  (5.27)  is  multiplied 
by  G  to  yield 

ei  =  GATn,  (5.30) 

where  the  vector  ei  is  a  normally  distributed,  zero  mean,  random  vector  with  covariance 
matrix 


Cei  =  £{GArnnTAGr} 

=  GCeoGT  (5.31) 

=  I. 

Thus  the  matrix  G  has  the  desired  effect  of  whitening,  a.k.a.  decorrelating,  the  error  vector 

ei. 

Returning  to  the  LP  representation  of  (5.23),  multiply  (5.23)  by  the  matrix  G  to 

yield 

ei  =  Gx0  +  GX0a0.  (5.32) 

The  error  power  of  (5.32),  | |ei 1 1|,  is  equal  to 

|  |ei  1 1 2  =  (Gx0  +  GX0a0)T(Gx0  +  GX0a0).  (5.33) 


It  can  be  shown  that  minimizing  the  expression  in  (5.33)  is  equivalent  to  minimizing  (5.24). 
If  G  is  assumed  not  to  be  a  function  of  aG,  the  weighted  least  squares  solution  to  minimize 
(5.33)  is 

a0  =  -(X^(GtG)X0)-1X^(GtG)x0.  (5.34) 

Inserting  the  value  of  G 1  G  and  Cn  =  a2I  from  (5.29)  into  the  above  equation  yields 

a0  =  -(X^(ATA)-1X0)-1X^(ArA)-1x0.  (5.35) 

Equation  (5.35)  correctly  accounts  for  the  coloring  of  the  equation  error  by  the  moving 
average  process. 
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G  is  a  function  of  aG.  However,  it  can  be  shown  that  sufficiently  close  to  aG  (5.35) 
is  a  contraction  mapping.  Thus,  when  (5.35)  is  iterated,  it  will  converge  to  a  fixed  point 
that  is,  in  view  of  (5.34)  and  (5.35),  close  to  the  minimum  of  the  error  given  by  (5.33) 
(and  thereby  (5.24))  -  see  [7,23-26].  Thus,  the  following  iterative  weighted  least  squares 
estimate  of  the  LP  coefficients  is  formed  from  (5.35) 

a£+1  =  —(X-o  (Afc  Afc)_1X0)~1XQ  (AjTAfc)-1x0,  (5.36) 

where  A/c  is  constructed  using  the  kth  iterate  of  the  LP  coefficients  aj,k\  and  a„fc+1')  of 
the  left  hand  side  of  (5.36)  is  the  k  +  1  iterate  of  the  LP  coefficients.  Equation  (5.32)  is 
minimized  when  a^+1^  =  a The  above  efficient  algorithm  is  referred  to  as  Iterative 
Generalized  Least  Squares  (IGLS).  A  good  initial  guess  for  the  values  of  a^  is  needed  to 
ensure  convergence  to  the  global  minimum,  especially  in  low  SNR.  Initializing  via  a  low- 
resolution  Periodogram  or  the  EPM  solution  is  normally  sufficient.  Note  that  initializing 
via  the  EPM  can  be  accomplished  by  setting  A1  A  =  I  for  the  first  iteration  of  (5.36). 
Simulation  experiments  have  established  10  iteration  steps  will  suffice. 

In  [7],  it  is  shown  that  minimizing  the  objective  function  via  (5.36)  is  equivalent 
to  obtaining  the  ML  estimate  of  the  LP  coefficients,  since  both  same  objective  function 
is  minimized.  By  the  invariance  property  of  ML  estimation,  since  there  is  a  one  to  one 
mapping  from  the  LP  coefficients  to  the  frequencies,  it  is  also  a  ML  estimate  of  the  frequen¬ 
cies  [11].  Although  derived  under  different  assumptions,  note  there  are  many  similarities 
between  the  IGLS  algorithm  and  the  Iterative  Quadratic  Maximum  Likelihood  (IQML) 
algorithm,  as  would  be  expected,  since  both  yield  ML  estimates  of  the  LP  coefficients. 
IQML  was  itself  shown  to  be  equivalent  to  the  Iterative  Pre-filtering  algorithm  of  Steiglitz 
and  McBride  [27].  For  brevity’s  sake,  those  similarities  are  not  discussed  here.  The  inter¬ 
ested  reader  should  refer  to  [14]  for  the  most  comprehensive  discussion  of  IQML  and  [8] 
for  a  comparison  of  IQML  and  IGLS.  The  author  notes  that  the  IGLS  algorithm  has  been 
successfully  applied  to  diverse  applications  [23,24,28]. 

5.2.3  Confidence  Intervals  for  IGLS  Estimates.  It  is  important  to  predict  the 
accuracy  of  the  frequencies’  estimates  for  EW  receivers,  especially  when  frequency  can 
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discriminate  different  threat  types  as  in  Fig.  1.3.  To  this  end,  confidence  intervals  must  be 
established.  Normally,  confidence  intervals  are  specified  in  terms  of  probabilities  of  error 
contained  within  k  standard  deviations 

P(—ka  <  e  <  ka )  =  1  —  p,  (5.37) 

where  a  is  the  estimate’s  standard  deviation  and  1  -p  is  the  probability  the  error  is  con¬ 
tained  in  the  confidence  interval.  Data  driven  confidence  intervals  for  IGLS  estimates  are 
developed  below. 

To  establish  confidence  intervals  for  the  LP  coefficients,  the  LP  coefficients’  variance 
must  be  determined.  This  covariance  can  be  determined  using  ML  theory  if  the  SNR  is 
above  threshold3.  Note  that  the  whitened  Linear  Prediction  equation  error  is 

ei  —  Gs0  =  GS0a0  +  GA1  n,  (5.38) 

where  S0  and  s0  are  the  signal  components  of  X0  and  x0  respectively.  Let  q  =  ei  —  Gs0. 
Since  (5.38)  is  a  linear  transformation  of  a  Gaussian  random  variable,  q  is  also  a  Gaussian 
random  variable 

q  =  fV(GS0a0, 1).  (5.39) 

Taking  the  log  of  the  pdf  of  q  (5.39)  and  then  the  gradient  with  respect  to  a0  yields 

=  S^GTq-S^GTGS0a0.  (5.40) 

US.Q 

Setting  the  right  hand  side  of  (5.40)  equal  to  zero  and  solving  for  aG  yields  the  Maximum 
Likelihood  estimate  for  the  LP  coefficients 


aD  ml  =  (So  GTGS0)_1Sj  Grq.  (5.41) 

The  covariance  of  the  LP  coefficients’  ML  estimation  error  (Note  that  this  estimator  is 
unbiased)  can  be  obtained  by  backtracking  to  Equation  (5.40)  and  factoring  its  right  hand 

3The  threshold  effect  is  discussed  in  Section  5.2.4. 
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side 


(5.42) 


=  (S^GtGSg)((Sq  GtGS0)_1Sq  G^q  —  ac). 

U3-o 

Note  the  form  of  (5.40)  matches  the  form  of  the  Cranrer-Rao  Lower  bound  theorem  given 
by  [10];  an  estimator  only  achieves  efficiency  if  the  partial  derivative  with  respect  to  the 
parameter  to  be  estimated  can  be  written  in  the  following  form 

^-W(l(«)-«),  (5.43) 

where  x  is  an  arbitrary  vector  of  interest,  J(0)  is  the  Fisher  information  matrix,  and  g(x)  is 
the  ML  estimator.  Thus,  the  Fisher  information  matrix  of  the  LP  coefficient  ML  estimate 
of  (5.41)  is  (S^G7GS0),  and  therefore  the  covariance  of  the  estimation  error,  provided 
ML  estimation  of  the  LP  coefficients,  is 

Ca„  =  (So  GtGS0)_1.  (5.44) 


Unfortunately,  (5.44)  features  SQ  which  is  not  available.  If  the  observation  matrix  is 
substituted  for  the  signal  matrix,  the  covariance  has  randomness  associated  with  it  because 
of  the  noise.  However,  if  the  SNR  is  high  enough  and  assuming  the  LP  coefficients  are 
sufficiently  close  to  the  true  value,  the  noise  can  be  neglected,  and  an  estimate  of  the  LP 
coefficient  estimation  error  covariance  matrix  is  [7] 


(Xo  GtGX0)_1 
°'2(Xo  (AtA)~1X0)”1. 


(5.45) 


The  CaD  diagonals  are  the  LP  coefficient’s  predicted  estimation  error  variance,  that  is 
Ca0(m,m)  =  £{(a0(m)  -  a0(m))2}. 

Since  above  the  threshold  the  LP  coefficient’s  estimation  error  covariance  matrix 
achieves  efficiency,  the  method  from  [11]  can  be  invoked  to  obtain  the  frequencies’  estima¬ 
tion  error  covariance  matrix 

J_1(f)  =  HTJ”1(a0)H,  (5.46) 
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where  the  matrix  H  is  defined  as  H(m,p)  =  ,  and  where  wp  is  the  mapping  from 

LP  coefficients  to  frequency.  Now,  the  unique  mapping  of  LP  coefficients  to  frequencies 
for  two  sinusoids  is4 


fP  =  wp(  a0)  =  i-cos 


-l 


— a(l)  +  (-l)*>+Va(l)2-4a(2)  +  8' 


(5.47) 


Let  gp  =  a(1)+('  1'> - VaA)  4a(2)+8^  ^ien  ugjng  ^he  above  equation,  the  following  deriva¬ 

tives  can  be  calculated  for  the  matrix  H  (values  of  derivatives  from  [7]) 


dwp(  aQ) 
daD{ 1) 


1 

^  \/f  -g'p 


(-i)Mi) 

\J  a(l)2  —  4a(2)  +  8 


dwP{&o)  =  1  2  (-!)p 

dao(2)  8tt^/i  -gj  _  xM1)2  -4«(2)  +8 


(5.48) 


(5.49) 


The  above  derivatives  coupled  with  (5.46)  yields 


Cf  =  HTCaoH, 


(5.50) 


where  Cf  is  the  frequency  estimation  error  covariance  the  values  of  the  matrix  H  are  defined 
in  (5.48)  and  (5.49).  The  estimated  error  variance  of  the  individual  frequency  estimates  lie 
along  the  diagonal  of  the  estimated  frequency  error  covariance,  C f(p,p)  =  £{(fp  ~  fp)2}- 

Because  of  the  non-linearity  of  the  transformation  in  (5.47),  it  is  very  difficult  to 
calculate  the  value  of  p  in  (5.37)  analytically.  Therefore,  Monte  Carlo  simulations  are 
employed  to  determine  the  percentage  of  estimates  within  one,  two  and  three  standard 
deviations  of  the  corresponding  predicted  variance  to  determine  confidence  intervals  above 
threshold. 


5.2.4  IGLS  Algorithm  Simulations.  MC  simulations  are  performed  to  validate 
the  IGLS  algorithm.  The  signal  (5.18),  with  two  real  sinusoids,  is  considered  in  all  of  the 

4  For  five  sinusoids  or  greater  a  closed  form  relationship  is  not  possible. 
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simulations.  In  the  simulations,  the  experimental  MSE  is  calculated  as 

1  M 

mse  =  m£(/-/*)2  (5-51) 

i—  1 

where  M  is  the  number  of  MC  trials  and  /)  is  the  ith  frequency  estimate.  The  negative  log  of 
the  MSE  is  plotted  to  highlight  good  performance.  The  SNR  is  defined  as  — 10 log10(2cr2). 
The  bias  is  calculated  as 

1  M 

bias  =  m  ~~  ■&)■  (5-52) 

i= 1 

The  emphasis  is  on  signals  with  low  SNR,  short  data  record,  and  close  frequencies  to 
demonstrate  the  robustness  of  the  IGLS  algorithm-based  parametric  receiver  in  difficult 
estimation  environments. 

Figure  5.1  contains  a  plot  of  MSE  vs.  SNR  for  a  Monte  Carlo  (MC)  simulation 
calculated  at  every  1  dB  using  the  EPM  estimation  algorithm  with  simulation  parameters: 
N=32,  [A!  =  1,  /i  =  0.227,  fa  =  f  ],  [A2  =  1,  h  =  0.207,  fa  =  f  ];  SNR  =  -101og10(2a2), 
M=1000.  Note  that  the  EPM  does  not  achieve  the  CRB. 

Figure  5.2  contains  a  plot  of  MSE  vs.  SNR  for  a  MC  simulation  calculated  at  every  0.5 
dB  using  the  IGLS  algorithm  with  the  same  parameters  as  Figure  5.1.  The  initial  estimate 
for  Figure  5.2  is  taken  from  the  extended  Prony  method  estimate  by  setting  AT A  =  I  in 
Equation  (5.36).  From  Figures  5.2(a)  and  5.2(b),  if  the  noise  threshold  is  defined  where 
the  experimental  MSE  is  within  2  dB  of  the  CR  bound,  then  the  noise  threshold  for 
/i  and  f2  is  10  dB  and  12  dB  respectively.  Below  the  noise  threshold,  the  algorithm’s 
performance  quickly  drops  well  below  the  CRB  because  the  noise  has  overwhelmed  the 
estimation  algorithm’s  ability  to  estimate  the  signal.  The  noise  threshold  is  inherent  to 
non-linear  estimation.  As  expected,  the  experimental  MSE  achieves  the  CRB  above  the 
noise  threshold  since  the  IGLS  algorithm  is  a  ML  estimator  of  frequency.  Below  the 
noise  threshold,  the  estimates  are  biased,  as  can  be  seen  in  Figures  5.2(c)  and  5.2(d)  thus 
comparison  to  the  CRB  is  not  appropriate.  However,  above  threshold  the  estimates  are 
unbiased,  and  CRB  comparison  is  appropriate.  Finally,  note  by  comparing  Fig.  5.3  to  Fig. 
5.1  that  the  IGLS  algorithm  significantly  outperforms  the  EPM  algorithm. 
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Figure  5.1:  Extended  Prony  Method  Estimation  Accuracy. 


Figure  5.3  contains  the  IGLS  algorithm  Monte  Carlo  simulation  results  for  the  same 
parameters  as  Figure  5.2,  except  the  initial  guess  is  now  assumed  to  be  obtained  from  the 
closest  32  point  FFT  sample  frequency  point  to  the  two  signals  frequencies.  This  FFT 
frequency  sample  point  is  used  for  both  initial  frequencies  guesses  and  is  converted  to  the 
initial  LP  coefficients  guess  by  using  Table  5.1.  From  Figures  5.3(a)  and  5.3(b),  if  the 
noise  threshold  is  defined  where  the  experimental  MSE  is  within  2  dB  of  the  CR  bound, 
then  the  noise  threshold  for  and  fi  is  4.5  dB  and  6  dB  respectively.  Thus,  a  better 
noise  threshold  is  obtained  by  using  more  information  about  the  frequencies  values.  As 
expected,  the  experimental  MSE  achieves  the  CR  bound  above  the  noise  threshold  since 
the  IGLS  algorithm  is  a  MLE  of  frequency.  Note  from  Figures  5.3(c)  and  5.3(d)  that  the 
estimate  becomes  unbiased  above  the  noise  threshold,  hence  comparison  to  the  CRB  above 
the  threshold  is  appropriate. 

Figure  5.4  contains  the  IGLS  algorithm  Monte  Carlo  simulation  results  for  the  same 
parameters  as  Figure  5.3  including  using  the  same  32  point  FFT  frequency  sample  point 
as  the  initial  guess,  except  the  number  of  sample  points  is  increased  to  M=128.  The  IGLS 
algorithm  noise  threshold  in  Figures  5.4(a)  and  5.4(b)  for  these  parameters  decreased 
dramatically  to  below  -5  dB.  Note  that  the  bias  in  Figures  5.4(c)  and  5.4(d)  is  negligible 
for  all  SNR’s  considered. 

Figure  5.5  contains  the  IGLS  algorithm  Monte  Carlo  simulation  results  for  the  same 
parameters  as  Figure  5.4,  except  the  initial  guess  is  from  the  EPM.  The  IGLS  algorithm 
noise  threshold  for  the  frequencies  in  Figures  5.5(a)  and  5.5(b)  is  7  dB  and  5  dB  respectively, 
much  higher  than  when  an  initial  guess  is  provided.  Thus  increasing  the  data  length  by  a 
factor  of  four  only  decreases  the  noise  threshold  by  5  dB  for  both  frequencies  when  using 
the  extended  Prony  method  to  initialize  the  estimate.  The  jagged  anomaly  below  the 
noise  threshold  contained  in  Figure  5.5(a)  is  partially  explained  by  zooming  in  around  the 
anomaly  in  both  the  bias  and  frequency  as  in  Figure  5.6  and  noticing  that  there  is  a  slight 
rise  in  bias  at  the  same  point  as  the  anomaly,  which  probably  means  the  IGLS  algorithm 
is  converging  to  a  local  minimum  more  often.  The  IGLS  MC  simulations  in  [8]  and  [7] 
exhibited  the  same  type  of  anomaly  below  the  threshold. 


5-16 
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Figure  5.3:  IGLS  Estimation  Accuracy  with  an  FFT  initial  guess  and  N=32 


(a)  /i  MSE 


(b)  f2  MSE 


Figure  5.4:  IGLS  Estimation  Accuracy  with  an  FFT  initial  estimate  and  N=128 


(a)  fi  MSE 


(b)  f2  MSE 


(c)  /i  Bias  (d)  f2  Bias 


Figure  5.5:  IGLS  Estimation  Accuracy  with  an  EPM  initial  estimate  and  N=128. 


x  Kf4 


(a)  /i  MSE  zoomed  (b)  /i  Bias  zoomed 

Figure  5.6:  Figure  5.5(a)  zoom  analysis:  Spike  area  of  Figure  5.5(a)  zoomed  in 

along  with  bias  of  Figure  5.5(b).  Note  the  dip  in  experimental  bias  which  coincides  with 
the  spike  reduction  in  experimental  MSE. 


5. 2. 4-1  IGLS  Confidence  Interval  Simulations.  Figure  5.7  contains  a  com¬ 
parison  of  the  average  predicted  LP  coefficient  estimation  error  variance  for  the  signal  of 
Fig.  5.1 

1  M 

Var pred{a0(m)}  =  —  ^  Qo(m,  m),  (5.53) 

1=1 

and  the  MC  experimentally  obtained  variance 

1  M 

Var  eXp{d0{m)}  =  -  -  —  -  ^(a0(m)j  -  mean  (a0(m)))2,  (5.54) 

1=1 

where  mean  is  calculated  as  ^  YllLi  ®o(m)-  Also  plotted  is  the  MSE  between  the  calculated 
variance  and  each  variance  estimate.  At  11  dB  for  both  a(l)  and  a( 2),  the  experimental 
variance  converges  to  the  calculated  variance  and  the  MSE  is  negligible.  At  this  point  the 
noise  power  can  be  neglected  and  (5.45)  becomes  a  good  estimate  of  the  variance. 

Fig.  5.8  establishes  confidence  intervals  for  frequency  estimates  via  Monte  Carlo 
simulation  for  the  signal  of  Fig.  5.1.  At  each  SNR  point,  the  number  of  frequency  estimates 
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(a)  a(l)  (b)  a(2) 

Figure  5.7:  Average  predicted  LP  estimation  error  variance,  MC  experimental  variance, 
and  MSE  between  MC  experimental  variance  and  predicted  LP  estimation  error  variance 
versus  SNR. 

within  K  a  is  calculated  by 

M 

TX\fi  -  f  I  <  k°i)  for  k  2>  3-  (5.55) 

7=1 

where  b,;  is  calculated  using  (5.50).  From  Fig.  5.8,  the  confidence  intervals  are  valid  for 
/i  and  f-2  greater  than  11  dB.  Rough  estimates  of  p  are  0.7,  0.95,  and  less  than  0.99 
for  K  =  1,  2,  3.  Thus,  good  data  driven  confidence  intervals  are  established  using  the 
measurements  and  knowledge  of  the  noise  variance. 

In  light  of  the  results  of  the  analysis  and  the  MC  simulations,  it  is  apparent  that 
the  IGLS  algorithm  yields  Maximum  Likelihood  estimates  of  frequencies  and  therefore, 
comparison  to  the  IDR-CRB  is  valid  above  the  noise  threshold. 

5.3  IGLS  algorithm-based  Parametric  Receiver  compared  to  the  IDR-CRB 

The  IGLS  algorithm-based  parametric  receiver,  which  is  shown  above  to  yield  Max¬ 
imum  Likelihood  estimates,  is  compared  to  the  IDR-CRB  by  employing  a  Monte  Carlo 
simulation.  The  simulation  signal  is  defined  for  these  simulations  as 

x{n)  =  Aicos(27r/it  +  cfii)  +  A2ftCos(27r/2i  +  </>2)  +  n(t).  n  =  0  . . .  N  —  1  (5.56) 
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(a)  A  (b)  A 

Figure  5.8:  Frequency  Estimate  Confidence  Intervals  for  K  =  1,  2,  and  3  u’s. 

where  A2b  is  the  minimum  value  of  A 2  for  the  ML  Estimator  to  still  attain  the  desired 
MSE  of  f2  as  determined  by  the  IDR-CRB.  The  simulation  analysis  follows  the  flowchart 
of  Figure  5.9.  First,  the  simulation  parameters  for  each  frequency  separation  of  interest  are 
input  into  the  Fig.  4.3  IDR-CRB  routine.  If  the  IDR-CRB  routine  finds  the  parameters  at 
any  frequency  separation  not  achievable,  the  parameters  are  not  considered  by  the  IGLS 
algorithm  Monte  Carlo  simulation  in  the  next  step.  If  the  parameters  are  viable,  the  IDR- 
CRB  routine  generates  the  A2  bound  amplitude,  A2b',  the  lowest  A2  amplitude  that  will 
still  theoretically  achieve  the  given  parameters.  The  IDR-CRB  is  plotted  for  valid  A2b 
values  when  this  routine  is  finished.  Also,  as  will  become  apparent  from  the  simulation 
results,  the  SNR  is  also  an  important  analysis  tool  and  is  defined  as 

SNRA2b  =  10 log  (H)  (5.57) 

The  SNRa2I>  is  also  plotted  after  the  IDR-CRB  routine.  The  achievable  parameters  along 
with  the  corresponding  A2b  values  are  sent  to  the  next  step,  which  is  a  Monte  Carlo 
simulation  of  the  signal  defined  in  Equation  (5.56).  At  each  frequency  separation  value 
A/,  the  corresponding  A2b  amplitude  is  assigned  along  with  the  signal  parameters  listed  in 
the  caption  and  1000  independent  IGLS  Monte  Carlo  trials  are  performed.  The  resulting 
experimental  MSE  and  bias  for  and  A  are  plotted.  For  a  successful  experimental 
confirmation  of  the  IDR-CRB,  the  f\  IGLS  estimate  experimental  MSE  should  be  less 
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Figure  5.9:  IDR-CRB  IGLS  Performance  Analysis 

Flowchart. 

than  the  desired  RMS  frequency,  /^cc,  and  the  f2  IGLS  estimate  experimental  MSE  should 
equal  the  desired  f^cc  for  all  frequency  separation  values. 

Figure  5.10  contains  the  MC  results  for  the  desired  RMS  frequency  estimation  ac¬ 
curacies  of  facc  =  and  facc  =  jjjjy  with  initialization  via  the  EPM  and  decreasing 
A /  evaluated  at  0.001  increments  for  the  parameters:  [A\  =  l,/i  =  0.207,  (j>i  =  0], 
[A2  =  A2b,f2  =  fi  -  A/,  02  =  A7/2  -  7r  -  Ar7r/i ] ;  SNR  =  101og10(^).  In  Figure  5.10(c), 
the  estimates  experimental  MSE  clearly  shows  the  IGLS  algorithm-based  parametric  re¬ 
ceiver  achieves  the  IDR-CRB  for  all  applicable  values  since  f\  MSE  is  below  the  desired 
MSE  at  all  points  and  the  f2  achieves  the  desired  MSE  within  1  dB  at  all  points.  In  Figure 
5.10(e),  the  experimental  bias  is  negligible,  thus  comparison  to  the  desired  MSE  is  appli¬ 
cable  for  both  frequency  estimates.  In  Figure  5.10(d),  the  estimates  MSE  does  not  achieve 
the  IDR-CRB  for  all  applicable  values.  Noting  that  the  SNR  of  A2b  of  Figure  5.10(b)  is 
well  below  the  SNR  of  A2b  of  Figure  5.10(a)  and  that  the  bias  is  no  longer  negligible  in 
5.10(f),  the  conjecture  is  the  noise  threshold  of  the  IGLS  algorithm  has  been  exceeded  for 
the  parameter  combination. 

Figure  5.11  contains  the  performance  analysis  results  for  the  same  parameters  as 
Figure  5.10  with  initialization  via  the  closest  frequency  sample  of  a  32  point  FFT.  Figure 
5.11(c)  matches  closely  to  Figure  5.10(c),  achieving  the  desired  MSE  within  1  dB  for  all 
applicable  points.  However,  the  f2  MSE  in  Figure  5.11(d)  achieves  the  desired  MSE  within 
4  dB  for  all  applicable  values,  which  means  the  estimate  is  right  at  the  FFT  initialized 
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IGLS  algorithm’s  noise  threshold  and  the  noise  threshold  is  lower  for  the  FFT  initialized 
IGLS  algorithm,  as  expected. 

Figure  5.12  contains  the  performance  analysis  results  with  the  measurement  samples 
increased  to  N  =  128,  the  desired  RMS  frequency  tightened  to  face  =  4^77  with  initial¬ 
izations  via  EPM  or  closest  32  point  FFT  frequency  sample  point.  The  IGLS  algorithm 
performs  as  expected,  once  again  validating  the  derived  IDR-CRB  in  the  FFT  initialization. 
The  noise  threshold  is  once  again  encountered  with  the  EPM  initialization5. 

Thus,  the  Instantaneous  Dynamic  Range  Cramer  Rao  Bound  algorithm  is  validated 
using  the  IGLS  algorithm.  From  the  experimental  results,  the  IDR-CRB  is  useful  for  tight 
bounds  on  parametric  receivers  performance.  However,  because  of  the  noise  threshold,  the 
IDR-CRB  does  not  provide  a  bound  that  can  be  used  for  parametric  receiver  analysis  for 
loose  frequency  estimation  requirements.  Thus,  in  the  next  section,  the  author  proposes 
a  method  to  determine  the  IDR  for  a  parametric  receiver  for  a  loose  frequency  estimate 
requirement. 

5.3.1  IGLS  Algorithm-based  Parametric  Receiver  IDR  for  Loose  Frequency  Esti¬ 
mates.  Loose  frequency  estimates  are  obtained  above  the  noise  threshold  for  a  paramet¬ 
ric  receiver.  When  the  amplitude  ratio  is  high,  so  that  the  low  amplitude  carrier  is  buried 
in  noise  and  the  measurement  is  basically  of  a  single  carrier,  overmodelling  becomes  an 
issue  that  causes  both  frequency  estimates  to  deteriorate  (the  CRB  does  not  consider  this 
issue)  as  can  be  seen  in  Fig.  5.10(d).  In  an  IGLS  algorithm-based  parametric  receiver,  for 
loose  frequency  estimates,  the  amplitude  ratio  at  the  point  where  the  estimate  begins  to 
deteriorate  can  be  considered  the  IDR  for  loose  frequency  estimate  requirements.  If  the 
lower  amplitude  signal  is  below  the  threshold,  it  is  better  to  reduce  the  model  order  to 
accurately  estimate  the  higher  amplitude  sinusoid.  Fig.  5.13  is  a  MC  simulation  with  the 
same  parameters  as  Fig.  5.2,  except  the  Sinusoid  1  SNR  is  fixed  at  23  dB  and  the  A 2  SNR 
is  increased  in  0.5  dB  steps  by  increasing  the  A 2  amplitude.  For  Sinusoid  2  SNR  below  7.5 
dB,  the  frequency  estimates  for  Sinusoid  1  are  almost  20  dB  below  efficiency,  even  though 

sThis  noise  threshold  will  be  encountered  even  by  functions  that  minimize  the  objective  functions  for 
the  two  frequency  estimates  directly.  It  is  inherent  to  non-linear  estimation. 
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Figure  5.10:  IGLS  performance  for  IDR-CRB  A2b  (EPM) 


(a)  IDR-CRB,  SNR  A2b,  FFT  init  (b)  IDR-CRB,  SNR  ,42b,  EPM  init 


(c)  fi  MSE,  FFT  init 


(d)  fi  MSE,  EPM  init 


(e)  fi  Bias,  FFT  init 


(f)  fi  Bias,  ,  EPM  init 


Figure  5.12:  IGLS  performance  for  IDR-CRB  A-2b  (M=128,  facc  =  7^). 


70 


Figure  5.13:  IDR  determination  for  two  sinusoids  in  white  noise. 

sinusoid  1  has  an  SNR  of  23  dB.  In  this  situation,  the  underlying  model  order  in  the  algo¬ 
rithm  should  be  reduced  to  efficiently  estimate  Sinusoid  1.  At  7.5  dB  both  sinusoids  are 
measured  efficiently,  thus  the  IDR  is  23  —  7.5  =  15.5  dB.  This  definition  and  simulation  of 
the  IGLS  parametric  algorithm’s  IDR  for  loose  frequency  estimates  along  with  the  IGLS 
algorithm  development  has  been  submitted  for  publication  [29]. 

5-4  Conclusion 

The  IGLS  frequency  estimation  algorithm  is  developed  and  shown  to  achieve  ML 
estimates.  The  IDR-CRB  is  then  experimentally  verified  using  the  IGLS  algorithm-based 
parametric  receiver  for  tight  frequency  estimate  requirements.  For  loose  requirements,  the 
IDR-CRB  is  shown  to  be  unachievable  due  to  the  noise  threshold  inherent  to  non-linear 
estimation.  Thus,  an  alternate  method  to  determine  the  IGLS  algorithm-based  parametric 
receiver  IDR  for  loose  estimation  accuracies  is  proposed  by  determining  the  threshold  point 
where  both  signals  are  measured  efficiently. 
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VI.  Conclusions 


6. 1  Introduction 

This  thesis  investigates  the  EW  receiver’s  IDR  for  the  following  IDR  definition:  IDR 
is  defined  as  the  maximum  signal  amplitude  ratio  for  a  given  frequency  estimation 
accuracy,  a  given  frequency  separation  and  a  given  SNR  [5].  The  measured  signal  is  as¬ 
sumed  to  consist  of  multiple  sinusoidal  tones  in  white  noise  that  fill  the  entire  measurement 
window,  i.e. 

x(n)  =  s(n)  +  w(n).  n  =  0...N  —  l  (6.1) 

Two  types  of  EW  receivers  are  considered;  a  DFT-Based  Intercept  Receiver  and  an  IGLS 
Algorithm-Based  Parametric  Receiver.  In  Chapter  III,  the  DFT-Based  Intercept  Receiver, 
using  the  novel  SLR  method,  is  evaluated  for  IDR.  The  results  provide  numerical  estimates 
for  57  dB  of  IDR.  The  resulting  analysis  has  a  direct  impact  on  digital  EW  receiver  analysis 
and  design.  In  Chapter  IV,  the  method  used  to  calculate  the  IDR-CRB  in  [1]  is  simplified, 
and  the  IDR-CRB  is  extended  to  real  signals.  In  Chapter  V,  the  novel  IGLS  frequency  es¬ 
timation  algorithm,  originally  researched  by  [8]  and  [7] ,  is  completely  developed  and  shown 
to  yield  ML  results.  The  IDR-CRB  validates  using  the  IGLS  Algorithm-Based  Parametric 
Receiver  for  tight  frequency  estimate  requirements.  For  loose  frequency  estimate  require¬ 
ments,  the  author  proposes  determining  the  parametric  receiver  IDR  based  on  when  both 
measurements  achieve  efficiency. 

6.2  Contributions 

In  Chapter  III,  the  main  contribution  is  a  solid  method  in  order  to  evaluate  a  DFT- 
based  intercept  receiver’s  IDR.  Also,  the  novel  SLR  method’s  IDR  is  evaluated  and  shown 
to  yield  the  IDR  in  Table  3.1.  This  IDR  is  independent  of  bin  spacing  for  signals  separated 
by  more  than  2  bins  and  does  not  have  the  associated  processing  gain  loss  and  widening 
of  the  main  beam  inherent  to  window  based  approaches. 

In  Chapter  IV,  an  improved  method  to  calculate  the  complex  IDR-CRB  is  introduced. 
The  old  method  in  [1]  relied  on  an  iterative  technique  that  is  difficult  to  implement.  The 
new  method  exploits  the  fact  that  the  Fisher  information  matrix  can  be  factored  as  shown 
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in  [6]  for  a  method  that  requires  no  iteration  and  is  quite  simple.  The  IDR-CRB  is 
extended  to  real  signals,  again  based  on  the  results  in  [6].  Plots  are  provided  to  evaluate 
the  results.  Finally,  it  is  firmly  established  that  a  receiver’s  frequency  estimates  should 
only  be  compared  to  the  IDR-CRB  if  the  estimates  are  unbiased,  which  is  generally  not 
the  case  in  intercept  receivers  but  is  true  for  most  parametric  receivers. 

In  Chapter  V,  the  novel  IGLS  frequency  estimation  algorithm,  developed  by  Dr. 
Pachter  and  researched  by  Zahirniak  [7]  and  Ingham  [8],  is  completely  developed  and 
shown  to  yield  ML  estimates.  The  author  feels  the  explanation  of  IGLS  provided  is  the 
most  concise  and  informative  to  date.  The  IDR-CRB  is  then  evaluated  using  the  IGLS 
Algorithm-Based  Parametric  Receiver.  For  tight  frequency  estimate  requirements,  the 
IDR-CRB  is  a  valid  bound  on  performance.  For  less  stringent  frequency  estimate  require¬ 
ments,  the  author  proposes  determining  the  IDR  based  on  when  both  estimates  achieve 
efficiency.  The  author’s  loose  definition  of  IDR  coupled  with  the  concise  and  informative 
description  of  IGLS  has  been  submitted  for  publication  [29] . 

6.3  Future  Work 

Areas  of  future  work  include 

•  Implement  the  SLR  method  results  using  a  lookup  table  and  compare  to  the  ideal 
results  in  Table  3.1. 

•  Use  the  method  of  Chapter  III  to  compare  straight  windowing  results  to  the  SLR 
method  results  for  IDR  in  Table  3.1. 

•  Consider  other  non-parametric  spectral/frequency  estimation  algorithms  for  the  in¬ 
tercept  receiver  such  as  the  minimum  variance  algorithm  described  in  [13]  and  com¬ 
pare  results  for  IDR  to  Table  3.1. 

•  Compare  other  algorithms  to  the  IDR-CRB  when  tight  frequency  estimates  are  re¬ 
quired. 

•  Apply  the  IGLS  algorithm  to  other  applications  such  as  a  discrete  frequency  rate 
estimator  in  a  software  radio. 
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•  Research  ways  to  improve  the  IGLS  noise  threshold  when  initialized  by  the  EPM, 
such  as  an  additional  contraint. 

6.4  Summary 

The  results  in  this  thesis,  whether  taken  collectively  or  individually,  represent  signif¬ 
icant  contributions  to  the  field  of  Electronic  Warfare.  The  IDR  analysis  performed  in  this 
thesis  standardizes  how  the  two  types  of  EW  receivers  IDR  should  be  evaluated.  These 
contributions  directly  impact  the  design  and  updates  for  digital  EW  receivers.  Thus,  these 
results  directly  and  positively  impact  USAF  operations  in  the  critical  EW  field. 
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Appendix  A.  Interference  in  an  EW  Receiver 

A .  1  Introduction 

As  with  any  system,  noise  interferes  with  the  EW  receiver  measurements.  Since  the 
noise  is  random,  the  noise  is  characterized  using  statistics  for  mathematical  analysis.  The 
statistical  noise  model  used  throughout  the  thesis  is  developed  below,  with  a  brief  review 
of  random  process  theory  preceding  the  noise  model  development. 

A. 2  Random  Process  Theory 

Random  signals  consist  of  an  ensemble  of  member  functions.  The  signal  measured  in 
an  experiment  may  only  be  a  single  member  of  a  large  ensemble  of  signals  with  a  certain 
probability  of  selection  associated  with  each  member  of  the  ensemble.  Thus,  instead  of 
a  deterministic  mathematical  description,  random  signals  are  specified  in  terms  of  their 
statistical  characteristics,  with  the  following  being  the  most  important: 

•  The  mean  of  a  random  signal  is  defined  as 

px{n)  =  £{x(n)}  .  (1.1) 

Thus,  the  mean  of  a  random  signal  at  time  n  is  the  expected  value  is  of  the  ensemble 
functions  at  time  n. 

•  The  autocorrelation  of  a  random  signal  is  defined  as 

Rxx(n,  m)  =  £  { x{n)x*(m )}  .  (1.2) 

where  Rxx(n,m )  is  an  autocorrelation  value.  A  signal  is  defined  to  be  Wide  Sense 
Stationary  (WSS)  if  the  autocorrelation  depends  only  on  the  time  difference,  not 
the  actual  time  (very  important)  and  the  mean  is  constant  over  time.  Note  that 
the  Fourier  transform  of  the  autocorrelation  function  of  a  WSS  process  is  the  Power 
Spectral  Density  (PSD)-  see  Appendix  B. 
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•  The  autocovariance  of  a  random  signal  is  defined  as 


Cxx{n ,  m)  =  Rxx(n,  m)  -  //x(n)//*(m)  (1.3) 

where  Cxx(n,m )  is  an  autocovariance  value. 

•  Note  that  the  autocorrelation  and  the  autocovariance  can  be  extended  to  cross¬ 
correlation  and  cross-covariance  between  two  separate  signals. 

The  above  concepts  are  applied  in  the  next  section  to  mathematically  model  thermal  noise. 

A. 3  Thermal  Noise 

In  1827,  the  British  Botanist  Robert  Brown  observed  that  small  pollen  grains  exhibit 
random  motion  in  water.  In  1905,  Albert  Einstein  (who  was  unaware  of  Brown’s  work) 
showed  that  small  particles  (around  10~4  cm)  move  randomly  due  to  the  constant  bom¬ 
bardment  from  the  molecules  of  the  medium  [30].  This  random  movement  is  now  called 
Brownian  motion.  In  1923,  Norbert  Wiener,  using  stochastic  theory,  derived  a  mathemat¬ 
ical  random  process  model  for  Brownian  Motion  called  the  Wiener  Process  [30]. 

The  Wiener  Process  is  used  extensively  to  model  the  random  motion  of  electrons 
in  electronic  devices  such  as  resistors,  capacitors,  inductors,  and  semiconductor  devices 
(which  make  up  the  RF  amplifier  of  the  EW  receiver).  These  random  fluctuations  in 
electron  density  interfere  with  the  information  bearing  signals  that  flow  through  these 
components.  This  interference  is  referred  to  as  thermal  noise.  As  the  number  of  devices 
increases,  the  collection  of  Wiener  Processes  become  a  zero-mean,  stationary  Gaussian 
Process  due  to  the  Central  Limit  Theorem  [30].  The  Gaussian  process  model  for  thermal 
noise  exhibits  a  PSD  that  is  flat  over  a  wide  range  of  frequencies,  and  is  referred  to  as 
Additive  White  Gaussian  Noise  (AWGN)  [30]. 

The  white  noise  PSD  is  [30] 

SwwU )  =  ^  Joules/Hz  (1.4) 
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where  N0  is  equal  to  kT  (where  k  is  boltzmann’s  constant,  and  T  is  the  temperature  in 
Kelvin  of  the  noise  source),  and  Swwif )  is  the  PSD  of  the  white  noise.  This  implies  that 
the  autocorrelation  function  which  is  the  IFT  of  the  PSD  is  a  delta  function, 

Rww{t)  =  cr25(t),  (1.5) 

where  a2  is  the  noise  variance.  From  (1.5),  each  realization  of  thermal  noise  in  time  is 
independent  from  the  next  realization.  For  the  general  case  of  colored  or  white  noise,  the 
distribution  of  real  noise  for  any  stationary  process  for  the  general  N- variate,  zero  mean 
Gaussian  case  is  given  as  [30] 

fw( w)  =  (2vr)^|R|“iexp  |-^wrR_1w|  ,  (1.6) 

where  the  matrix  R  is  the  covariance  matrix  of  the  noise  (for  white  noise  R  =  cr2I,  where  I 
is  the  identity  matrix)  and  w  is  the  AWGN  vector.  Thermal  noise  is  assumed  throughout 
the  thesis. 
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Appendix  B.  Power  Spectral  Density  and  the  Periodogram 


B.  1  Introduction 


Although  the  DFS  is  used  for  the  frequency  estimation  algorithm  in  Section  3.3,  the 
Periodogram  -  an  estimate  of  a  WSS  process’s  PSD  1  -  is  still  an  important  tool 
for  frequency  estimation,  and,  because  the  Periodogram  is  a  power  statistic,  has  a  natural 
relation  to  the  direct  Maximum  Likelihood  estimate  of  multiple  frequencies.  Thus,  in 
Section  B.2,  the  deterministic  PSD  is  developed.  In  Section  B.3,  the  PSD  of  a  WSS 
random  signal  is  introduced  and  then  estimated  with  the  Periodogram. 


B.2  Deterministic  Power  Spectral  Density 

If  an  arbitrary  signal  voltage  is  referenced  to  a  one  ohm  resistor,  then  the  energy  of 
the  signal  is  defined  as  [19] 

E  =  lim  f  \x(t)\2dt,  (2-1) 

T^oo  J_T 

and  the  power  of  the  signal  is  defined  as 

T 

p  ~  1™  f  [l  \x(t)\2dt.  (2.2) 

It  is  a  well-known  relationship  [16, 19]  that  energy  in  the  frequency  domain  equals  the 
energy  in  the  time  domain 

/OO  POO 

\x(t)\2dt=  /  \X(f)\2df.  (2.3) 

-oo  J — OO 

Equation  (2.3)  is  known  as  Parseval’s  theorem.  A  similar  property  applies  to  the  power 
relationship  between  the  time  and  frequency  representation  of  a  signal  [19] 

1  /■?  1  r°° 

p=  lim  -  /  \x(t)\2dt  =  lim  -  /  \X(f)\2df.  (2.4) 

OO  T  J  =T  T—kx  T  J _00 

1  Since  the  signal’s  contained  in  the  measurement  signal  are  deterministic,  the  measurement  is  actually 
a  mixed  process,  this  fact  is  normally  ignored  in  frequency /spectral  estimation. 


B-l 


For  periodic  signals  the  limiting  operation  can  be  ignored  in  the  power  equations  of  above, 
and  Ta  can  be  substituted  for  T,  where  T0  is  the  period  of  the  signals.  A  similar  relationship 
that  was  made  in  (2.3)  can  be  applied  in  the  discrete  time  and  frequency  domains 

TV— 1  TV-1 

P=^2\x(n)\2  =  -^2\X(k)\2.  (2.5) 

n=0  k= 0 


Equation  (2.4)  can  be  rewritten  as 

1  /'  -  /“OO 

p  =  lim  -  /  \x(t)\2dt  =  /  Sx(f)df.  (2.6) 

T^oo  1  J=T  .7-00 

where  Sx(f)  is  the  PSD  of  a  signal  defined  as 

Sx(f)=lunJ-^pl.  (2.7) 

The  discrete  PSD  is  defined  as 


Sx[k]  =  lim  hx(k)\2. 

M— xx)  iv 


(2.8) 


where  N  is  the  number  of  discrete  frequency  sample  points.  The  PSD  represents  the  power 
per  unit  Hertz  of  a  signal.  An  interesting  result  is  the  equivalent  operation  in  the  time 
domain  to  the  PSD 


1 

F~1{SX}=  lim  -  /  \X(f)\2e^adf 

1  — XX  1 


'  —  OO 

poo 


lim  -  /  X(f)X(f)*eP^adf 

^oo  1  J ^ 

T 

X(f)(  f  ^  x{t)e-j2nftdt)*ej2nfadf 

lim  -  /  x*(t)  /  X(f)ej2nPt+a)dfdt 

"^OO  T  Jt  J 


=  lim  — 

T— >oo  T 


'  — OO 

r*oo 


’  — oo 
.  T 


„  T 

1  f* 

=  lim  —  /  x  (t)x(t  +  a)dt 

OO  T  J_T 


=  Rf(a), 


(2.9) 
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where  Rf{oc)  is  the  autocorrelation  function  and  has  a  central  role  in  random  process 
theory. 

The  PSD  can  be  employed  in  digital  EW  systems  for  the  detection/estimation  of 
sinusoids.  Since  the  period  of  the  sinusoids  is  unknown,  the  measurement  time  (window 
time  length)  is  substituted  for  the  period.  Thus,  the  PSD  of  the  signal  in  (3.3)  is 

sx(f)  =  WS-.  (2.10) 

r 

To  make  the  task  easier,  rewrite  X(f)  as 

X(f)  =  AB  +  CD  (2.11) 

where  A,B,C,  and  D  are  defined  as  follows, 

A  —  TJ^1  p-jftr/lT+f ) 

2 

B  =  sine (t(/  +  -  sine (r(/  -  fi))ej‘t>1 

n  —  l^lp-j(7rhT+  f ) 

2 

D  =  sine (t(/  +  /2))e-^2  -  sine (r(/  -  f2))eJ<t>2 
Then  the  PSD  of  x{t )  is 

SJf)  =  ~(\AB\2  +  A*  B*  CD  +  ABC*D*  +  \CD\2).  (2.16) 

T 

Although  the  above  formula  is  for  the  continuous  frequency  case,  it  will  match  closely  to 
the  discrete  PSD  at  the  frequency  sample  points  with  slight  differences  due  to  the  aliasing 
of  the  sine  function  side-lobes.  The  cross  terms  of  Equation  (2.16)  increase  the  frequency 
estimates  bias  for  closely  spaced  sinusoids.  The  deterministic  PSD  plot  results  mirror  the 
previous  DFT  magnitude  plots  in  Chapter  III  for  most  results  (since  those  results  were 
reported  in  dB),  except  for  the  jj  scaling  term. 


(2.12) 

(2.13) 

(2.14) 

(2.15) 


B-3 


B.3  Periodogram 

For  WSS  random  signals,  the  PSD  is  defined  as  [13] 

Sxx(f)  =  jjm  g||Xi^/)|2|  •  (2.17) 

Ignoring  the  expected  value  operator  and  using  the  available  data  yields  the  PSD  estimate 

Sxx(f)  =  IXrWI\  (2.18) 

T 

where  Sxx(f)  is  called  the  Periodogram.  The  Periodogram  is  an  inconsistent  spectral 
estimator,  i.e.  the  variance  does  not  decrease  as  the  number  of  measurements  increases. 
The  Periodogram  is  often  employed  for  frequency  detection/estimation.  However,  because 
the  interpolation  algorithms  in  Section  3.3  rely  on  phase  information,  the  Periodogram  is 
not  used  in  this  thesis. 
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